> M <- 10000
> ctr <- 0
> n1 <- 20
> n2 <- 5    # unequal sample sizes
> tstar <- qt(0.975, df = n1+n2-2)
> 
>   for (mc in 1:M){
+     Y1 <- rexp(n1, rate = 1/sqrt(2)) - sqrt(2)
+     Y2 <- runif(n2,-0.5,0.5)       # rt(n2, df=4)    
+     sp <- sqrt(var(Y1)*(n1-1) + var(Y2)*(n2-1)) / sqrt(n1+n2-2)
+     half <- tstar*sp*sqrt(1/n1 + 1/n2)
+     if (abs(mean(Y2) - mean(Y1)) > half){
+         ctr <- ctr + 1 
+     }  
+   }
> 
>  print( paste("failure rate of supposedly 95% C.I.", ctr/M) ) 
[1] "failure rate of supposedly 95% C.I. 0.0102"
> 
> 
>  X <- read.csv("U:/public_html/588/Class4/cloud.csv", header =TRUE)
> 
> dim(X)
[1] 52  2
> X[1:5,]
  RAINFALL TREATMENT
1   1202.6  UNSEEDED
2    830.1  UNSEEDED
3    372.4  UNSEEDED
4    345.5  UNSEEDED
5    321.2  UNSEEDED
> boxplot(X[,1] ~ X[,2])
> hist(X[1,])
Error in hist.default(X[1, ]) : 'x' must be numeric
> hist(X[,1])
> Z <- log(X[,1])
> hist(X)
Error in hist.default(X) : 'x' must be numeric
> hist(Z)
> boxplot(Z ~ X[,2])
> t.test(Z ~ X[,2])

        Welch Two Sample t-test

data:  Z by X[, 2] 
t = 2.5444, df = 49.966, p-value = 0.01408
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 0.2408498 2.0467125 
sample estimates:
  mean in group SEEDED mean in group UNSEEDED 
              5.134187               3.990406 

> exp(c(0.2408, 2.047))
[1] 1.272267 7.744632
> qqnorm(Z)
>