# sim_CI is adapted from the function conf.int() written by Yihui Xie, see:
# http://finzi.psych.upenn.edu/R/library/animation/html/conf.int.html
#
# sim_CI calculates confidence intervals (CIs) with confidence level "alpha"
# (default: .95) of "samp" samples (default: 50) of samples size "n" (default:
# 10) drawn from a population with a theoretical mean "mean" (default: 0) and a
# standard deviation "sd" (default: 1) based on the t- or normal distribution
# (default: t-distribution with t.dist=T) and the probability of a one- or two-
# sided test (default: two-sided with one.sided=F).
#
# The CIs will be plotted (default: with plot=T), the color of CIs not including
# the actual population mean will be red. The speed of drawing the CIs can be
# adjusted (default: .25 secs per CI with interval=.25).
#
# The seed of the random generator can be set to any value to create replicable
# data using "seed=" (default: no seed).
#
# The resulting object (obtainable with "object=sim_CI()") will contain the
# alpha level ($alpha), the sample size ($n), the empirical means ($mean) and
# standard deviations ($sd) of all samples, the lower and upper limits of the
# CIs ($CI), the distribution used to calculate the CIs ($distr) and the
# coverage rate ($CR), i.e. the rate of confidence intervals including the
# actual population mean.
#
# Examples see:
# http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Software/CI_demo.r

sim_CI = function(mean=0,sd=1,n=10,alpha=.95,samp=50,t.dist=T,one.sided=F,
                  plot=T,interval=0.25,seed=NULL)
{
  if (length(seed)!=0) set.seed(seed)
  d = replicate(samp,rnorm(n,mean,sd))
  m = apply(d,2,mean)
  s = apply(d,2,sd)
  level = ifelse(one.sided,alpha,1-(1-alpha)/2)
  test = ifelse(t.dist,qt(level,n-1),qnorm(level))
  CI.L = m - test*s/sqrt(n)
  CI.U = m + test*s/sqrt(n)
  rg = c(mean-max(abs(c(mean-CI.L,mean-CI.U))),
         mean+max(abs(c(mean-CI.L,mean-CI.U))))
  cvr = (CI.L < mean) & (CI.U > mean)
  if (plot==T) {
    xax = pretty(1:samp)
    if (t.dist==T) {
      if (one.sided==T) {
        mtit = expression("CI = "~bar(x)%+-%t[alpha]~"*"~hat(sigma)/sqrt(n)~"")
      }
      else {
        mtit = expression("CI = "~bar(x)%+-%t[alpha/2]~"*"~hat(sigma)/sqrt(n)~"")
      }
    }
    else {
      if (one.sided==T) {
        mtit = expression("CI = "~bar(x)%+-%z[alpha]~"*"~hat(sigma)/sqrt(n)~"")
      }
      else {
        mtit = expression("CI = "~bar(x)%+-%z[alpha/2]~"*"~hat(sigma)/sqrt(n)~"")
      }
    }
    plot(1:samp,ylim=rg,type="n",xaxt = "n",
         xlab=paste("Sample (n=",n,", µ=",mean,", sigma=",sd,")",sep=""),
         ylab=paste(100*alpha,"%-CI",sep=""),main = mtit)
    abline(h=mean,lty=2)
    for (i in 1:samp) {
      axis(1, xax[xax <= i])
      if (i==1) {Sys.sleep(2)}
      arrows(i, CI.L[i], i, CI.U[i], length=par("din")[1]/samp * 0.5,
             angle=90, code=3, col=c("red","gray")[cvr[i] + 1])
      points(i, m[i])
      erg=legend("topright", legend=table(factor(cvr[1:i],levels=c('FALSE','TRUE'))),
                 fill=c("red","gray"), bty="n")
      Sys.sleep(interval)
      if (i != samp) {
      rect(erg$rect$left,erg$rect$top-0.015*erg$rect$h,
           erg$rect$left+erg$rect$w-0.015*erg$rect$w,erg$rect$top-erg$rect$h,
           border="white",col="white")
      }
    }
  }
  CI=cbind(CI.L,CI.U)
  if (one.sided==T) {
    colnames(CI) = paste(round(c((1-alpha), 1-(1-alpha))*100, 2), "%")
  }
  else {
    colnames(CI) = paste(round(c((1-alpha)/2, 1-(1-alpha)/2)*100, 2), "%")
  }
  rownames(CI) = 1:samp
  distr = ifelse(t.dist,"t-distribution","normal-distribution")
  invisible(list(alpha=alpha,n=n,mean=m,sd=s,CI=CI,distr=distr,CR=mean(cvr)))
}