0

Below I have my code, which determines the number of times the mean of a population falls within confidence intervals of samples taken from the mean. Basically, hoping to prove 95% confidence intervals work.

rp<-function(x,s,n){ #x-population data, s-number of samples taken from 
population, n-size of samples
  m<-mean(x) 
  ci.mat=NULL
  tot=0
  for(i in 1:s){
   cix<-t.test(sample(x,n))$conf.int #obtain confidence intervals of 1000 
samples of x
    if(cix[1]<m & m<cix[2]){tot<-tot+1} #total number of confidence intervals containing pop mean
    ci.mat<-rbind(ci.mat,cbind(cix[1],cix[2]))
  }
  par(mfrow=c(2,1))
  hist(ci.mat[,1],main=paste("Lower Limits for Sample Confidence 
Intervals"),xlab="Lower Limit")
  hist(ci.mat[,2], main=paste("Upper Limits for Sample Confidence 
Intervals"),xlab="Upper Limit")
  return(data.frame(mean(x),tot/s))
}

I am hoping to add the population mean to my histograms, so I can show where the confidence intervals did not include the mean. So wherever the mean lies on the histogram with the lower limits, any values to the right of it would be part of a confidence interval that did not include the mean. I have no experience with modifying plots in R, so I don't even know if this is possible. Thanks for any help!

Lost
  • 331
  • 1
  • 12
  • So what would the desired output look like? It would help if you provided a [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input data so we can run the code. Maybe sketch what you would like the output to look like so we can have a better understanding of what you are after. – MrFlick Oct 16 '17 at 19:04
  • I just used simulated normal data, so an example of running the code would be: rp(rnorm(10000,1,1),1000,100) What I am looking for in the histograms is to show which lower limits were too high to include the population mean, and which upper limits were too low to include the pop mean. – Lost Oct 16 '17 at 19:07
  • Try calling `abline(v = m)` after each `hist(.)`. Also, since you change `par(mfrow)`, I would **start** the function with the following two lines of code: `op <- par(mfrow=c(2,1)); on.exit(par(op))`. (Two lines, here separated by a semi-colon but two lines.) – Rui Barradas Oct 16 '17 at 19:32
  • Thank you, Rui Barradas! abline was what I needed! – Lost Oct 16 '17 at 20:59

0 Answers0