# How to understand confidence level? K <- 100 alpha <- 0.1 n <- 10 mu <- 15 sigma <- 2 CI_L <- CI_U <- rep(NA, K) for (k in 1:K){ samples <- rnorm(n, mu, sigma) xbar <- mean(samples) s <- sd(samples) qtile <- qt(1-alpha/2, n - 1) CI_L[k] <- xbar - qtile * s / sqrt(n) CI_U[k] <- xbar + qtile * s / sqrt(n) } plot(1:K, type = 'n', ylim = c(min(CI_L), max(CI_U)), xlab = '', ylab = paste0((1-alpha)*100, '%CI')) segments(1:K, CI_L, 1:K, CI_U, col = 2, lwd = 2) abline(h = mu, lty = 'dashed') sum(CI_L <= 15 & CI_U >= 15)/K ### Confidence interval CI for mu when sigma is unkown muCI <- function(x, alpha = 0.05, sigma0 = 1, unknown = TRUE) { n <- length(x) mu.val <- mean(x) sd.val <- sd(x) cri.val <- qt(1 - alpha/2, n - 1) if (unknown == FALSE) { sd.val <- sigma0 cri.val <- qnorm(1 - alpha/2, 0, 1) } CI <- c(mu.val - sd.val * cri.val/sqrt(n), mu.val + sd.val * cri.val/sqrt(n)) list(CI = CI, conflevel = 1 - alpha) } ## eg x <- c(4.68, 4.85, 4.32, 4.85, 4.61, 5.02, 5.2, 4.6, 4.58, 4.72, 4.38, 4.7) muCI(x, alpha = 0.05) muCI(x, alpha = 0.01) # CI for sigam2 when sigma is unknown sigma2CI <- function(x, alpha = 0.05) { n <- length(x) sd2.val <- var(x) CI <- c((n - 1) * sd2.val/qchisq(1 - alpha/2, n - 1), (n - 1) * sd2.val/qchisq(alpha/2, n - 1)) list(CI = CI, conflevel = 1 - alpha) } ## eg x <- c(45.3, 45.4, 45.1, 45.3, 45.5, 45.7, 45.4, 45.3, 45.6) sigma2CI(x) ### CI for parameter for two normal distributions musigma2CI <- function(x, y, alpha = 0.05, sigma10 = 1, sigma20 = 1, unknown = TRUE) { m <- length(x) n <- length(y) mu.x <- mean(x) mu.y <- mean(y) s2.x <- var(x) s2.y <- var(y) sd2.w <- ((m - 1) * s2.x + (n - 1) * s2.y)/(m + n - 2) sd.val <- sqrt(sd2.w) * sqrt(1/m + 1/n) cri.mu <- qt(1 - alpha/2, m + n - 2) cri.sigma2 <- qf(c(1 - alpha/2, alpha/2), m - 1, n - 1) if (unknown == FALSE) { sd.val <- sqrt(sigma10^2/m + sigma20^2/n) cri.mu <- qnorm(1 - alpha/2, 0, 1) } CI.mu <- c(mu.x - mu.y - sd.val * cri.mu, mu.x - mu.y + sd.val * cri.mu) CI.sigma2 <- s2.x/s2.y * 1/cri.sigma2 list(conflevel = 1 - alpha, CI.mu = CI.mu, CI.sigma2 = CI.sigma2) } ## eg x <- c(15, 14.8, 15.2, 15.4, 14.9, 15.1, 15.2, 14.8) y <- c(15.2, 15, 14.8, 15.1, 14.6, 14.8, 15.1, 14.5, 15) musigma2CI(x, y, alpha = 0.1, sigma10 = 0.18, sigma20 = 0.24, unknown = F) musigma2CI(x, y, alpha = 0.1, unknown = T)