setwd("D:/Courses/Statistical Computing/Data/")
setwd("E:/Courses/Statistical Computing/Data/")

rain.data <- read.csv("Seeding_Expt.csv",
header=T,na.strings = "-9")

attach(rain.data)

rain.data


# Set up two rainfall vectors, one for "no effect" and one for "10% increase".

rain.no <- TelDan
rain.10 <- rain.no-(seed==0)*rain.no/11


#  Summary when seeding has no effect.

tapply(rain.no,seed,mean)

T1 <- mean(rain.no[seed==1])-mean(rain.no[seed==0])
T2 <- mean(rain.no[seed==1])/mean(rain.no[seed==0])

c(T1,T2)


#  Summary when seeding has 10% effect.

tapply(rain.10,seed,mean)

T1 <- mean(rain.10[seed==1])-mean(rain.10[seed==0])
T2 <- mean(rain.10[seed==1])/mean(rain.10[seed==0])

c(T1,T2)


boot.2sample <- function(B,y0,y1){
  n0 <- length(y0)
  n1 <- length(y1)
  T1.boot <- NULL
  T2.boot <- NULL
  for (b in 1:B){
    y0.b <- sample(y0,n0,replace=TRUE)
    y1.b <- sample(y1,n1,replace=TRUE)
    T1.boot <- c(T1.boot,mean(y1.b)-mean(y0.b))
    T2.boot <- c(T2.boot,mean(y1.b)/mean(y0.b))
  }
  print("Summary stats for T1")
  print(c(mean(T1.boot),sd(T1.boot)))
  print(c(mean(T2.boot),sd(T2.boot)))
  par(mfrow=c(1,2))
  hist(T1.boot)
  hist(T2.boot)
  return(c(mean(T1.boot),sd(T1.boot),mean(T2.boot),sd(T2.boot)))
}

output <- boot.2sample(400,rain.no[seed==0],rain.no[seed==1])

output <- boot.2sample(400,rain.10[seed==0],rain.10[seed==1])



boot.2sample.ci <- function(B,y0,y1,alpha){
  n0 <- length(y0)
  n1 <- length(y1)
  T1.boot <- NULL
  T2.boot <- NULL
  for (b in 1:B){
    y0.b <- sample(y0,n0,replace=TRUE)
    y1.b <- sample(y1,n1,replace=TRUE)
    T1.boot <- c(T1.boot,mean(y1.b)-mean(y0.b))
    T2.boot <- c(T2.boot,mean(y1.b)/mean(y0.b))
  }
  print("Normal 1-alpha CI using T1")
  print(mean(T1.boot)-qnorm(1-alpha/2)*sd(T1.boot))
  print(mean(T1.boot)+qnorm(1-alpha/2)*sd(T1.boot))
  print("Normal 1-alpha CI using T2")
  print(mean(T2.boot)-qnorm(1-alpha/2)*sd(T2.boot))
  print(mean(T2.boot)+qnorm(1-alpha/2)*sd(T2.boot))
  print("Percentile 1-alpha CI using T1")
  print(quantile(T1.boot,alpha/2))
  print(quantile(T1.boot,1-alpha/2))
  print("Percentile 1-alpha CI using T2")
  print(quantile(T2.boot,alpha/2))
  print(quantile(T2.boot,1-alpha/2))
  return(c(mean(T1.boot),sd(T1.boot),mean(T2.boot),sd(T2.boot)))
}

output <- boot.2sample.ci(400,rain.no[seed==0],rain.no[seed==1],0.10)

output <- boot.2sample.ci(400,rain.10[seed==0],rain.10[seed==1],0.10)