# R/confIntDemo.R In CarletonStats: Functions for Statistics Classes at Carleton College

#### Documented in confIntDemo

```#' Confidence Interval Demonstration
#'
#' Draw many random samples and compute confidence interval.  How many
#' intervals capture the true mean?
#'
#' This simulation will draw 100 random samples from a given population
#' distribution and compute the correpsonding confidence intervals. The 100
#' intervals will be drawn with an indication of the ones that missed the true
#' mean. A histogram of the population will also be created.
#'
#' @param distr distribution of the population to be sampled. Options include
#' \code{"normal"}, \code{"exponential"}, \code{"uniform"} and \code{"binary"}
#' (partial match allowed).
#' @param size sample size
#' @param conf.level confidence level.
#' @return The command invisibly returns the fraction of intervals that capture
#' the true mean.
#' @author Laura Chihara
#' @keywords confidence interval
#' @examples
#'
#' confIntDemo()
#'
#' confIntDemo(distr = "exponential", size = 40)
#'
#'
#' @importFrom stats rbinom rnorm rexp runif qt qnorm qt sd
#' @export
confIntDemo <-
function(distr = "normal", size = 20, conf.level = .95)
{

if (size > 100000) stop("Choose smaller sample size")
if (size < 2) stop("Choose larger sample size")

alpha <- 1 - conf.level
conf.text <- deparse(substitute(conf.level))
title.text<-paste("100 confidence intervals, confidence level", conf.text,sep=" ")

popdist <- pmatch(distr, c("binary", "normal", "exponential", "uniform"), nomatch = NA)
if (is.na(distr))  stop("Distribution must be one of \"binary\", \"normal\", \"exponential\" or \"uniform\"")

pop <- switch(popdist,
binary = rbinom(1000000, 1, .6),
normal = rnorm(1000000, 10, 2),
exponential = rexp(100000, rate = 10),
uniform = runif(1000000, 0, 20))

#plot population distribution
if (popdist != 1)
{

hist(pop, xlab = "", main = "")
title("Population")
}

mu <- mean(pop)

index <- logical(100)
count <- 0
temp <- matrix(0,  nrow = 100, ncol = 2)
for (i in 1:100)
{
x <- sample(pop, size, replace = F)
xbar <- mean(x)
if (popdist == "binary")
{se <- sqrt(xbar*(1 - xbar)/size)
lower <- xbar-abs(qnorm(alpha/2))*se
upper <- xbar+abs(qnorm(alpha/2))*se
}
else
{
se <- sd(x)/sqrt(size)
lower <- xbar-abs(qt(alpha/2, size - 1))*se
upper <- xbar+abs(qt(alpha/2, size - 1))*se
}
temp[i,1] <- lower
temp[i,2] <- upper

if (lower < mu && mu < upper)
{count <- count + 1
index[i] <- T

}
}

xmin <- min(temp[,1])
xmax <- max(temp[,2])

#Rstudio doesn't have a dev.new() command
if (Sys.getenv("RSTUDIO_USER_IDENTITY") == "")
grDevices::dev.new()

plot(seq(xmin,xmax,length = 101), 0:100, type = "n", xlab = "", ylab = "")
abline(v = mu)

for (j in 1:100)
{
lines(temp[j,], c(j,j), col = ifelse(index[j], 1, 2), lwd = ifelse(index[j], 1, 2))
}

if (popdist == 1)
{ mtext("True proportion", side = 3, at = mu)
sub.text <- paste(count, " intervals cover true proportion", sep = " ")
title(main = title.text, sub = sub.text)

}
else
{mtext("True mu", side = 3, at = mu)
sub.text <- paste(count, " intervals cover true mean", sep = " ")
title(main = title.text, sub = sub.text)
}

invisible(count/100)

}
```

## Try the CarletonStats package in your browser

Any scripts or data that you put into this service are public.

CarletonStats documentation built on Aug. 10, 2018, 1:14 a.m.