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)

jack.2sample <- function(y,group){
  n <- length(y)
  T1.jack <- NULL
  T2.jack <- NULL
  for (i in 1:n){
    if (i==1){
      y1 <- y[2:n]
      gr1 <- group[2:n]
    }
    if (i==n){
      y1 <- y[1:(n-1)]
      gr1 <- group[1:(n-1)]
    }
    if ((i>1)&(i<n)){
      y1 <- c(y[1:(i-1)],y[(i+1):n])
      gr1 <- c(group[1:(i-1)],group[(i+1):n])
    }   
    T1.jack <- c(T1.jack,mean(y1[gr1==1])-mean(y1[gr1==0]))
    T2.jack <- c(T2.jack,mean(y1[gr1==1])/mean(y1[gr1==0]))
    i <- i+1
  }
  print("Summary stats for T1")
  print(c(mean(T1.jack),((n-1)*sqrt(1/n))*sd(T1.jack)))
  print(c(mean(T2.jack),((n-1)*sqrt(1/n))*sd(T2.jack)))
  return(T1.jack,T2.jack)
}


jack.2sample(rain.no,seed)