set.seed(123)
## true values
mu <- 0.0
sd <- 1.0

n <- 10 # sample size
nbatches <- 1000 # number of samples


## Functions that computes the probability that a reference value trueValue
## falls into confidence intervals computed by sampling from a normal with
## mean center and standard deviation stdev.
getProbConf <- function(center, stdev, trueValue, alpha, n, nbatches){
  # center   : mean of the normal from which is sampled
  # stdev    : standard deviation of the normal
  # trueValue: mean of the distribution from which sample was drawn
  # alpha    : significance level
  # n        : sample size
  # nbatches : number of samples that will be drawn
  q <- qnorm(p=1-alpha/2)
  tmp <- t(replicate(n=nbatches, expr={
    rn <- center+stdev*rnorm(n)
    ci <- mean(rn) + c(-1,1)*q*stdev/sqrt(n)
    findInterval(x=trueValue, vec=ci)
  }))
  return(mean(tmp == 1))
}


## Probability that mu falls into confidence intervals computed from samples
## from a normal N(mu,sd). Should be ~0.95
getProbConf(center=mu, stdev=sd, trueValue=mu,
            alpha=0.05, n=n, nbatches=nbatches)


## We draw some observation from the true distribution
x <- mu+sqrt(sd)*rnorm(n)
x_bar <- mean(x)
sd_hat <- sd(x)


## Probability that x_bar falls into confidence intervals computed from samples
## from a normal N(x_bar,sd_hat). Should be ~0.95
getProbConf(center=x_bar, stdev=sd_hat, trueValue=x_bar,
            alpha=0.05, n=n, nbatches=nbatches)




## Probability that mu falls into confidence intervals computed from samples
## from a normal N(x_bar,sd_hat).
getProbConf(center=x_bar, stdev=sd_hat, trueValue=mu,
            alpha=0.05, n=n, nbatches=nbatches)


## Do the same for nrep samples from N(mu, sd). Takes 10 seconds...
nrep <- 1000
probs <- replicate(nrep, {
  x <- mu+sqrt(sd)*rnorm(n)
  x_bar <- mean(x)
  sd_hat <- sd(x)
  getProbConf(center=x_bar, stdev=sd_hat, trueValue=mu,
              alpha=0.05, n=n, nbatches=nbatches)
})

plot(density(probs))
## mean
mean(probs)
## mode
d <- density(probs)
d$x[which.max(d$y)]

## On average, the true mean (mu) from which a sample with mean x_bar and
## standard deviation sd_hat was generated will not lie with a probability of
## 0.95 in a 0.95-confidence interval constructed from a sample
## from N(x_bar, sd_hat)