#-----------------------------------------------------#
#---   Classical and Bayesian inference (AMS132)   ---#
#---   Class 17                                    ---#
#-----------------------------------------------------#


# confidence intervals
mu <- 5
s2 <- 2.5
s <- sqrt(s2)
n <- 26
gamma <- 0.9

c <- qt((1+gamma)/2, df=n-1)

#  sample one:
set.seed(1)
x <- rnorm(n, mean=mu, sd=s)

xbar <- mean(x) # sum(x)/n
sprime <- sd(x) # sqrt( sum((x - xbar)^2)/(n-1) )

confint <- c(xbar - c*sprime/sqrt(n), xbar + c*sprime/sqrt(n))

#  sample two:
set.seed(3)
x <- rnorm(n, mean=mu, sd=s)

xbar <- mean(x) # sum(x)/n
sprime <- sd(x) # sqrt( sum((x - xbar)^2)/(n-1) )

confint <- c(xbar - c*sprime/sqrt(n), xbar + c*sprime/sqrt(n))



# several samples:
reps <- 100
confints <- matrix(0, ncol=2, nrow=reps)
counts <- rep(0,reps)
for(i in 1:reps){
  set.seed(4*i)
  x <- rnorm(n, mean=mu, sd=s)
  
  xbar <- mean(x) # sum(x)/n
  sprime <- sd(x) # sqrt( sum((x - xbar)^2)/(n-1) )
  
  confints[i, ] <- c(xbar - c*sprime/sqrt(n), xbar + c*sprime/sqrt(n))
  
  if(mu > confints[i, 1] & mu < confints[i, 2]){ counts[i] <- 1  }
}
mean(counts)

# plot:
plot(0,0, xlim=c(3.5, 7), ylim=c(0,100))
for(i in 1:reps){
  segments(confints[i,1], i, confints[i,2], i, lwd=2); readline()
  points(mu, i, pch=20, col="red"); readline()
}