R : Copyright 2002, The R Development Core Team
Version 1.5.1  (2002-06-17)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> invisible(options(echo = TRUE))
> foo <- read.table(url("http://www.stat.umn.edu/geyer/5601/hwdata/t6-6.txt"),
+     header = TRUE)
> attach(foo)
> names(foo)
[1] "information" "number"     
> 
> dat <- number
> grp <- as.ordered(information)
> mu0 <- mean(dat)
> sig0 <- sd(dat)
> n <- length(dat)
> ss.null <- sum((dat - mu0)^2)
> 
> xbar <- sapply(split(dat, grp), mean)
> nbar <- sapply(split(dat, grp), length)
> k <- length(xbar)
> 
> pava <- function(x, w) {
+     if (any(w <= 0))
+         stop("weights must be positive")
+     if (length(x) != length(w))
+         stop("arguments not same length")
+     n <- length(x)
+     design <- diag(n)
+     repeat {
+         out <- lm(x ~ design + 0, weights = w)
+         mu <- coefficients(out)
+         dmu <- diff(mu)
+         if (all(dmu >= 0)) break
+         j <- min(seq(along = dmu)[dmu < 0])
+         design[ , j] <- design[ , j] + design[ , j + 1]
+         design <- design[ , - (j + 1), drop = FALSE]
+     }
+     return(as.numeric(design %*% mu))
+ }
>  
> test.stat <- function(x, w) {
+     mu <- pava(x, w)
+     mu0 <- sum(w * x) / sum(w)
+     ss.alt <- sum(w * (x - mu)^2)
+     ss.null <- sum(w * (x - mu0)^2)
+     return(ss.null - ss.alt)
+ }
> 
> print(tstat <- test.stat(xbar, nbar))
[1] 52.77778
> 
> invisible(runif(1)); print(.Random.seed)
[1]         1 409518761 117368298
>  
> nsim <- 1e4 - 1
> tsim <- double(nsim)
> for (i in 1:nsim) {
+     xbarsim <- rnorm(k, mu0, sig0 / sqrt(nbar))
+     tsim[i] <- test.stat(xbarsim, nbar)
+ }
> phat <- mean(tsim >= tstat)
> (nsim * phat + 1) / (nsim + 1)
[1] 0.0358
> nsim / (nsim + 1) * sqrt(phat * (1 - phat) / nsim)
[1] 0.001855408
> 
> proc.time()
[1] 365.31   6.38 373.68   0.00   0.00
>