###################################################
### 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: CInormal
###################################################
n <- 50
alpha <- .05
UCL <- replicate(1000, expr = {
    x <- rnorm(n, mean = 0, sd = 3)
    (n-1) * var(x) / qchisq(alpha, df = n-1)
    } )


###################################################
### chunk number 3: CInormal2
###################################################
sum(UCL > 9)


###################################################
### chunk number 4: CInormal3
###################################################
mean(UCL > 9)


###################################################
### chunk number 5: CIlognormal
###################################################
sigma <- 1.125
UCL.logn <- replicate(1000, expr = {
    x <- rlnorm(n, meanlog = 0, sdlog = 1.125)
    (n-1) * var(x) / qchisq(alpha, df = n-1)
    } )
v <- exp(sigma^2)*(exp(sigma^2) - 1)
v
print(p <- mean(UCL.logn > v))


###################################################
### chunk number 6: 
###################################################
    ucl <- UCL.logn[1:100]
    plot(ucl, 1:100)
    segments(0, 1:100, ucl, 1:100, col=4)
    abline(v=9)


###################################################
### chunk number 7: 
###################################################
 normix <- function(n, M, S, p) {
     u <- runif(n)
     dst <- ifelse(u <= p, 1, 2)
     mu <- M[dst]
     sigma <- S[dst]
     x <- rnorm(n, mean=mu, sd=sigma)
     }


###################################################
### chunk number 8: CInormix
###################################################
UCL.mix <- replicate(1000, expr = {
    x <- normix(n, c(0,3), c(3,3), p=.5)
    (n-1) * var(x) / qchisq(alpha, df = n-1)
    } )
print(p <- mean(UCL.mix > 9))


###################################################
### chunk number 9: testH0
###################################################
n <- 50
alpha <- .05
sigma0 <- 1.125
x <- rlnorm(n, 0, sigma0)
t0 <- (n-1)*var(x)/sigma0^2
tr <- replicate(1000, expr = {
    x <- rlnorm(n, 0, sigma0)
    (n-1)*var(x)/sigma0^2
    } )
t0
quantile(tr, .95, type=1)


###################################################
### chunk number 10: 
###################################################
mean(tr >= t0)


###################################################
### chunk number 11: testH0hist
###################################################
plot(density(tr), xlim=c(0,10000),
   main="Ex. 4a: Replicates of T")
points(sum(t0>=tr), 0, pch=17, cex=2)


###################################################
### chunk number 12: testH0
###################################################
n <- 50
alpha <- .05
sigma0 <- 1.125; sigma1 <- 1.5
x <- rlnorm(n, 0, sigma1)
t0 <- (n-1)*var(x)/sigma0^2   # under H0
tr <- replicate(1000, expr = {
    x <- rlnorm(n, 0, sigma0)
    (n-1)*var(x)/sigma0^2
    } )
t0
quantile(tr, .95, type=1)


###################################################
### chunk number 13: 
###################################################
mean(tr >= t0)


###################################################
### chunk number 14: testH0hist2
###################################################
plot(density(tr), main="Ex. 4b: Replicates of T")
points(sum(t0>=tr), 0, pch=17, cex=2)