# R/CIsim.R In BSDA: Basic Statistics and Data Analysis

#### Documented in CIsim

```#' Confidence Interval Simulation Program
#'
#' This program simulates random samples from which it constructs confidence
#' intervals for one of the parameters mean (Mu), variance (Sigma), or
#' proportion of successes (Pi).
#'
#' Default is to construct confidence intervals for the population mean.
#' Simulated confidence intervals for the population variance or population
#' proportion of successes are possible by selecting the appropriate value in
#' the type argument.
#'
#' @param samples the number of samples desired.
#' @param n the size of each sample.
#' @param mu if constructing confidence intervals for the population mean or
#' the population variance, mu is the population mean (i.e., type is one of
#' either \code{"Mean"}, or \code{"Var"}). If constructing confidence intervals
#' for the poulation proportion of successes, the value entered for mu
#' represents the population proportion of successes \code{(Pi)}, and as such,
#' must be a number between 0 and 1.
#' @param sigma the population standard deviation. \code{sigma} is not required
#' if confidence intervals are of type \code{"Pi"}.
#' @param conf.level confidence level for the graphed confidence intervals,
#' restricted to lie between zero and one.
#' @param type character string, one of \code{"Mean"}, \code{"Var"} or
#' \code{"Pi"}, or just the initial letter of each, indicating the type of
#' confidence interval simulation to perform.
#' @return Graph depicts simulated confidence intervals. The number of
#' confidence intervals that do not contain the parameter of interest are
#' counted and reported in the commands window.
#' @author Alan T. Arnholt
#' @keywords distribution
#' @examples
#'
#' CIsim(100, 30, 100, 10)
#'     # Simulates 100 samples of size 30 from
#'     # a normal distribution with mean 100
#'     # and standard deviation 10.  From the
#'     # 100 simulated samples, 95% confidence
#'     # intervals for the Mean are constructed
#'     # and depicted in the graph.
#'
#' CIsim(100, 30, 100, 10, type="Var")
#'     # Simulates 100 samples of size 30 from
#'     # a normal distribution with mean 100
#'     # and standard deviation 10.  From the
#'     # 100 simulated samples, 95% confidence
#'     # intervals for the variance are constructed
#'     # and depicted in the graph.
#'
#' CIsim(100, 50, .5, type="Pi", conf.level=.90)
#'     # Simulates 100 samples of size 50 from
#'     # a binomial distribution where the population
#'     # proportion of successes is 0.5.  From the
#'     # 100 simulated samples, 90% confidence
#'     # intervals for Pi are constructed
#'     # and depicted in the graph.
#'
#' @export CIsim
CIsim <-
function(samples=100, n=30, mu=0, sigma=1, conf.level = 0.95, type ="Mean")
{
Aorange <- "#FF4C0C"
alpha <-1-conf.level
CL<-conf.level*100
N <-samples
choices <- c("Mean", "Var", "Pi")
alt <- pmatch(type, choices)
type <- choices[alt]
if (length(type) > 1 || is.na(type))
stop("alternative must be one \"Mean\", \"Var\", \"Pi\"")
if (type == "Pi" && (mu <=0 |mu >= 1))
stop("Value for Pi (mu) must be between 0 and 1.")
if (N <= 0 ||  n <= 0)
stop("Number of random CIs (samples) and sample size (n) must both be at least 1")
if (!missing(conf.level) && (length(conf.level) != 1 || !is.finite(conf.level) || conf.level <= 0 || conf.level >= 1))
stop("'conf.level' must be a single number between 0 and 1")
if (sigma <= 0 &&  (type=="Var" || type=="Mean") )
stop("Variance must be a positive value")

if (type == "Mean")
{
junk <- rnorm(N*n, mu, sigma)
jmat <- matrix(junk, N, n)
xbar <- apply(jmat, 1, mean)
ll <- xbar - qnorm(1 - alpha/2)*sigma/sqrt(n)
ul <- xbar + qnorm(1 - alpha/2)*sigma/sqrt(n)
notin <- sum((ll > mu) + (ul < mu))
percentage <- round((notin/N) * 100,2)
plot(ll, type = "n", ylim = c(min(ll), max(ul)), xlab = " ", ylab = " ")
title(sub=bquote(paste("Note: ",.(percentage),"% of the random confidence intervals do not contain ", mu ,"=", .(mu))))
title(main=bquote(paste(.(N), " random ", .(CL), "% confidence intervals where ", mu,  " = ", .(mu) )))
for(i in 1:N)
{
low<-ll[i];
high<-ul[i];
if(low < mu & high > mu)
{
segments(i,low,i,high)
}
else if(low > mu & high > mu )
{
segments(i,low,i,high, col=Aorange, lwd=5)
}
else
{
}
}
abline(h = mu)
cat(percentage,"% of the random confidence intervals do not contain Mu =", mu,".", "\n")
}
else if (type == "Var")
{
junk <- rnorm(N*n, mu, sigma)
jmat <- matrix(junk, N, n)
s2 <- apply(jmat, 1, var)
ll <- ((n - 1)*s2)/qchisq(1 - alpha/2, (n - 1))
ul <- ((n - 1)*s2)/qchisq(alpha/2, (n -1))
variance <- sigma^2
notin <- sum((ll > variance) + (ul < variance))
percentage <- round((notin/samples) * 100,2)
plot(ll, type = "n", ylim = c(min(ll), max(ul)), xlab = " ", ylab = " " )
title(sub=bquote(paste("Note: ",.(percentage),"% of the random confidence intervals do not contain ", sigma^2 ,"=", .(variance))))
title(main=bquote(paste(.(N), " random ", .(CL), "% confidence intervals where ", sigma^2,  " = ", .(variance) )))
for(i in 1:N)
{
low<-ll[i]
high<-ul[i]
if(low < variance & high > variance)
{
segments(i,low,i,high)
}
else if( low > variance & high > variance )
{
segments(i,low,i,high, col=Aorange, lwd=5)
}
else
{
}
}
abline(h = variance)
cat(percentage,"% of the random confidence intervals do not contain Var =", sigma^2,".", "\n")
}
else if (type == "Pi")
{
X <- rbinom(samples, n, mu)
p <- X/n
ll <- p - qnorm(1 - alpha/2)*sqrt((p * (1 - p))/n)
ul <- p + qnorm(1 - alpha/2)*sqrt((p * (1 - p))/n)
notin <- sum((ll > mu) + (ul < mu) )
percentage <- round((notin/samples)*100,2)
plot(ll, type = "n", ylim = c(min(ll), max(ul)), xlab = " ", ylab = " " )
title(sub=bquote(paste("Note: ",.(percentage),"% of the random confidence intervals do not contain ",pi,"=",.(mu))))
title(main=bquote(paste(.(N), " random ", .(CL), "% confidence intervals where ", pi,  "=", .(mu) )))
for(i in 1:N)
{
low<-ll[i]
high<-ul[i]
if( low < mu & high > mu)
{
segments(i,low,i,high)
}
else if( low > mu & high > mu )
{
segments(i,low,i,high, col=Aorange, lwd=5)
}
else
{