Confidence Interval Simulation

#1
Hello all,

I am trying to write a confidence interval simulation in R. The idea is to simulate several samples from a known distribution, calculate the mean and 95% CI for the mean, and then plot all confidence intervals in such a way that every confidence interval catching the "unknown parameter" will be in one color and others in different color. I also wish to calculate the percentage of CI's "catching" the parameter.

I have stared doing it:

Code:
n = 30
numOfSamples = 10
mu = 0
std = 2


d = data.frame()

for (sampleID in 1:numOfSamples)
{
  x = rnorm(n,mu,std)
  s = rep(sampleID,n)
  temp = cbind(x,s)
  d = rbind(d,temp)
}


for (sampleID in 1:numOfSamples)
{
  d_temp = subset(d,s==sampleID)
  avg = mean(d_temp$x)
  ci = t.test(d_temp$x, data = d_temp)$conf.int
}
What I need your help with is:

1. Saving the means and CI's of each sample in a dataset.
2. Creating the plot, which should look like this:
Plot
or something similar. I prefer ggplot2 for it to be nicer, but any plot will do....

I believe that my code is nearly ready, just saving and plotting... :D

Thank you in advance !
 

Dason

Ambassador to the humans
#2
This is some really old code I wrote a long time ago that does something similar. I didn't really read your question in too much depth but hopefully this helps?

Code:
CiSim2<-function(n,alpha,m,...){
  #n - sample size
  #alpha - alpha level for confidence regions
  #m - number of intervals to simulate
  #... - additional arguements passed to plot, such as xlim, ylim, etc...
  low <-1:m
  high <- 1:m
  captured <- rep(TRUE,m)
  c1 <- qnorm(1-alpha/2)/sqrt(n)
  c12 <- c1^2
  
  for(j in 1:m){
    x <- rnorm(n)
    xbar <- mean(x)
    s <- sd(x)
    captured[j] <- xbar^2 >= c12*s^2
    low[j] <- xbar - c1*s
    high[j] <- xbar + c1*s
  }
  
  plot(c(low,high),type="n",xlim=c(1,100),xlab="",ylab="",pch=19,...)
  abline(h = 0, lty = 1,lwd=3)
  points(high, col = 1, pch = 20)
  points(low, col = 1, pch = 20)
  
  cols <- c("black", "red")
  for(i in 1:100){
    lines(c(i,i), c(low[i],high[i]), lty = captured[i]+1, col = cols[captured[i]+1], pch = 19)
  }
  title(paste(100*(1-alpha),"% with n = ",n))
}

CiSim2(30, .05, 100)

par(mfrow=c(2,1))
#CiSim2(30,.20,100,ylim=c(-0.75,0.75))
CiSim2(30,.05,100,ylim=c(-.9,.9))
#CiSim2(120,.2,100,ylim=c(-0.75,0.75))
CiSim2(120,.05,100,ylim=c(-.9,.9))