###################################################
### chunk number 1: startup
###################################################
    rm(list=ls())
    set.seed(20090804)
    options(width=55,digits=4,continue=" ")
    par(ask=FALSE, mfrow=c(1,1))


###################################################
### chunk number 2: gof
###################################################
R <- 999
n <- 50
p.vals <- replicate(R, expr={
    x <- rnorm(n)
    shapiro.test(x)$p.value })
mean(p.vals <= .10)       #should be apprx. .10


###################################################
### chunk number 3: gof1
###################################################
R <- 999
n <- 50
p.vals <- replicate(R, expr={
    x <- rgamma(n, shape=4, rate=1)
    shapiro.test(x)$p.value })
mean(p.vals <= .10)       #power estimate


###################################################
### chunk number 4: power.sw
###################################################
 R <- 2000; n <- 50
 k <- seq(5, 100, 5)   #for large k, apprx. normal
 K <- length(k)
 pwr <- numeric(K)
 for (i in 1:K) {
    p <- replicate(R, expr={
        x <- rgamma(n, shape=k[i], rate=1)
        shapiro.test(x)$p.value})
    pwr[i] <- mean(p <= .10)
 }


###################################################
### chunk number 5: swpower2
###################################################
library(Hmisc)
se <- sqrt(pwr*(1-pwr)/R)
plot(k, pwr, type="l", ylim=c(0,1), ylab="power",
     main="SW power against Gamma(k,1), n=50 at 10%")
points(k, pwr, cex=.5)
errbar(k, pwr, pwr+2*se, pwr-2*se, add=TRUE)
detach(package:Hmisc)


###################################################
### chunk number 6: sleep
###################################################
 extra <- sleep$extra
 group <- sleep$group
 extra
 group


###################################################
### chunk number 7: t-statistic
###################################################
  t0 <- t.test(extra ~ group)$statistic   #theta-hat
  t0


###################################################
### chunk number 8: t-replicates
###################################################
  R <- 999
  t.star <- replicate(R, expr={
    g <- sample(group)  #a permutation of group labels
    t.test(extra ~ g)$statistic })


###################################################
### chunk number 9: ASL
###################################################
(sum(abs(t.star) >= abs(t0)) + 1) / (R+1)  #ASL or p-value


###################################################
### chunk number 10: permdist
###################################################
hist(t.star, main="", prob=TRUE, nclass="scott")
points(t0, 0, pch=17, cex=2)


###################################################
### chunk number 11: t-test
###################################################
t.test(extra ~ group)