cont=scan("ozonecontrol.csv")
treat=scan("ozonetreat.csv")
boxplot(list(cont, treat), names = c("control", "treatment"))
stripchart(list(cont, treat), group.names=c("control" "treatment"))
qqplot(cont, treat); grid()
t.test(cont,treat)

#permutation test
n = length(cont)
m=length(treat)
N=m+n
x = c(cont,treat)
dd =function(x,m,N){
	id =sample(1:N,m)
	dd = mean(x[id])-mean(x[-id])}
	
tst = replicate(10000,dd(x,m,N))
hist(tst,25)
dobs = abs(mean(cont)-mean(treat))
sum(abs(tst) > dobs)/10000  #p-value

#mann-whitney (rank sum) test
	
w = wilcox.test(cont,treat)
w$statistic/(length(cont)*length(treat)) #pi-hat
w$p.value

#bootstrap pi-hat
pihatstar = replicate(1000,wilcox.test(sample(cont,n,replace=T),
   + sample(treat,m,replace=T))$statistic)/(m*n)
hist(pihatstar,20)
sd(pihatstar)

#bootstrap ci for difference of means

dboot = replicate(10000, mean(sample(cont,n,replace=T))- mean(sample(treat,m,replace=T)))
mean(dboot)
delta = mean(cont)-mean(treat)
delta
sd(dboot)
sqrt(var(cont)/n + var(treat)/m)
lo = quantile(dboot-delta, .025)
hi = quantile(dboot-delta, .975)
 c(delta - hi, delta - lo) #the bootstrap confidence interval -- see p 284-85

#compare to ci based on t distribution
t.test(cont,treat)$conf.int