# Confidence Interval Simulation

#### Yankel

##### Member
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...

Thank you in advance !

#### Dason

##### Ambassador to the humans
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))

#### Dason

##### Ambassador to the humans
That code isn't entirely correct. And it pains me to look at it. I'll try doing a better rewrite later.