#Tutorial-3

#Finite Population CLT
yp<-islands####The areas in thousands of square miles of the landmasses which
#exceed 10,000 square miles.##################
mu<-mean(yp)
sig2<-var(yp)
N<-length(yp)
n<-5
s<-sample(1:N,n)
ys<-yp[s]

### Usual CLT holds for SRS with replacement
simmeanwr<-function(yp,n,b) {
  N<-length(yp)
  ybar<-c()
  for (i in 1:b) {
    s<-sample(1:N,n,replace=T)
    ys<-yp[s]
    ybar[i]<-mean(ys)
  }
  var<-((N-1)/N)*(sig2/n)
  out<-list(ybar,var)
  return(out)
}

b<-10000
par(mfrow=c(2,3))
for (n in c(5,10,30,50,100,200)) {
  sim<-simmeanwr(yp,n,b)
  ybardist<-(sim[[1]]-mu)/sqrt(sim[[2]])
  title<-toString(n)
  hist(ybardist,freq=F,main=title)
  z<-seq(min(ybardist),max(ybardist),length=100)
  lines(z,dnorm(z),col="red")
}


#### Finite population CLT for SRS without replacement
simmean<-function(yp,n,b) {
  N<-length(yp)
  sig2<-var(yp)
  ybar<-c()
  for (i in 1:b) {
    s<-sample(1:N,n)
    ys<-yp[s]
    ybar[i]<-mean(ys)
  }
  var<-((N-n)/N)*(sig2/n)
  out<-list(ybar,var)
  return(out)
}

b<-10000
for (n in c(5,10,15,24,35,45)) {
  sim<-simmean(yp,n,b)
  ybardist<-(sim[[1]]-mu)/sqrt(sim[[2]])
  hist(ybardist,freq=F,ylim=c(0,1.2))
  z<-seq(min(ybardist),max(ybardist),length=100)
  lines(z,dnorm(z),col="red")
}



#We are going to estimate the exact coverage probability for the approximate (1-alpha)% CIs through simulation
setwd("C:/Users/Shirin/Desktop/STUFF/STAT-410")
catsData<-read.table(file="catsM.csv",header=T,sep = ",") 
y.p<-catsData$Hwt#population
mu<-mean(y.p) #population mean
N<-length(y.p) #population size
alpha<-.05 #so that (1-alpha) is .95
n<-10 #sample size
b<-10000 #number of iterations in simulation
counter<-0
lower<-c()
upper<-c()
for (i in 1:b) {
  y.s<-sample(y.p,n)
  ybar<-mean(y.s)
  var.hat.ybar<-((N-n)/N)*var(y.s)/n
  delta<-qnorm((1-alpha/2))*sqrt(var.hat.ybar)
  lower[i]<-ybar-delta
  upper[i]<-ybar+delta
  if (mu>lower[i] & mu<upper[i]) counter<-counter+1 #check if the CI covers the populayion mean
}
cover.prob<-counter/b #estimated coverage probability
#Do the same with a larger n; what do you expect the results to be?

cov.p<-function(y.p,n){
  mu<-mean(y.p)
  counter<-0
  lower<-c()
  upper<-c()
  for (i in 1:b) {
    y.s<-sample(y.p,n)
    ybar<-mean(y.s)
    var.hat.ybar<-((N-n)/N)*var(y.p)/n
    delta<-qnorm((1-alpha/2))*sqrt(var.hat.ybar)
    lower[i]<-ybar-delta
    upper[i]<-ybar+delta
    if (mu>lower[i] & mu<upper[i]) counter<-counter+1 #check if the CI covers the populayion mean
  }
  cover.prob<-counter/b
}

cp<-c()
j<-0
for (n in seq(10,90,length=20)) {
  j<-j+1
  cp[j]<-cov.p(y.p,n)
}
par(mfrow=c(1,1))
plot(cp,ylim=c(.92,.98),pch=16)
abline(h=.95,col="red")