Lab Homework 5 - Confidence Interval and Hypothesis Tests

1. Let’s say a random varible \(X\) follows a normal distribution of \(N(0, 1)\).

(a) Draw 5 samples of \(X\) and use \(\bar{x}\) to infer the true population mean using a confidence interval estimate with confidence level of 95%. Do a simulation of 10,000 times and find the probability of the estimated interval actually contains \(\mu = 0\).

# CI for normal mean
ci_mean <- function(x, conf = 0.95) {
  n <- length(x)
  m <- mean(x)
  s <- sd(x)
  error <- qt(1 - (1 - conf)/2, df = n - 1) * s / sqrt(n)
  c(lower = m - error, upper = m + error)
}

n_sim <- 10000
n_samp <- 5
mu <- 0

ci_contains_mu <- replicate(n_sim, {
  x <- rnorm(n_samp, mean = mu, sd = 1)
  ci <- ci_mean(x, conf = 0.95)
  mu >= ci[1] && mu <= ci[2]
})
prob_ci_contains_mu <- mean(ci_contains_mu)
cat("Probability CI contains mu:", prob_ci_contains_mu, "\n")

Probability CI contains mu: 0.9516 

(b) Do a hypothesis test for \(X\) with \(H_0: \mu = 0\) and \(H_a: \mu \neq 0\) with 10,000 simulations and \(\alpha = 0.05\). Find the probability of rejecting \(H_0\) by comparing the p-value and \(\alpha\).

reject_H0_twosided <- replicate(n_sim, {
  x <- rnorm(n_samp, mean = mu, sd = 1)
  t.test(x, mu = mu, alternative = "two.sided")$p.value < 0.05
})
prob_reject_H0_twosided <- mean(reject_H0_twosided)
cat("Probability of rejecting H0 (two-sided):", prob_reject_H0_twosided, "\n")

Probability of rejecting H0 (two-sided): 0.0509 

(c) Redo (b) with \(H_a: \mu > 0\).

reject_H0_greater <- replicate(n_sim, {
  x <- rnorm(n_samp, mean = mu, sd = 1)
  t.test(x, mu = mu, alternative = "greater")$p.value < 0.05
})
prob_reject_H0_greater <- mean(reject_H0_greater)
cat("Probability of rejecting H0 (greater):", prob_reject_H0_greater, "\n")

Probability of rejecting H0 (greater): 0.0494 

(d) Increase the number of samples to 100 and repeat (a)-(c).

n_samp_large <- 100

ci_contains_mu_large <- replicate(n_sim, {
  x <- rnorm(n_samp_large, mean = mu, sd = 1)
  ci <- ci_mean(x, conf = 0.95)
  mu >= ci[1] && mu <= ci[2]
})

prob_ci_contains_mu_large <- mean(ci_contains_mu_large)
cat("Probability CI contains mu (n=100):", prob_ci_contains_mu_large, "\n")

Probability CI contains mu (n=100): 0.952 

reject_H0_twosided_large <- replicate(n_sim, {
  x <- rnorm(n_samp_large, mean = mu, sd = 1)
  t.test(x, mu = mu, alternative = "two.sided")$p.value < 0.05
})
prob_reject_H0_twosided_large <- mean(reject_H0_twosided_large)
cat("Probability of rejecting H0 (two-sided, n=100):", prob_reject_H0_twosided_large, "\n")

Probability of rejecting H0 (two-sided, n=100): 0.0532 

reject_H0_greater_large <- replicate(n_sim, {
  x <- rnorm(n_samp_large, mean = mu, sd = 1)
  t.test(x, mu = mu, alternative = "greater")$p.value < 0.05
})
prob_reject_H0_greater_large <- mean(reject_H0_greater_large)
cat("Probability of rejecting H0 (greater, n=100):", prob_reject_H0_greater_large, "\n")

Probability of rejecting H0 (greater, n=100): 0.0489 

(e) How does this experiment verify the theory that we learn? Summarize your findings properly.

Answer: - For a 95% confidence interval, the empirical coverage is close to 0.95, confirming the theoretical expectation. - The probability of rejecting \(H_0\) is approximately 0.05, matching the significance level (\(\alpha = 0.05\)). - Increasing the sample size improves the accuracy of the confidence interval coverage, as predicted by the Central Limit Theorem.

2. Repeat all steps in problem 1 with \(X\) follows a uniform distribution \(\text{Unif}(0,1)\). How does the distribution of \(X\) affect the results?

# For Unif(0,1), true mean is 0.5
mu_unif <- 0.5

ci_contains_mu_unif <- replicate(n_sim, {
  x <- runif(n_samp, min = 0, max = 1)
  ci <- ci_mean(x, conf = 0.95)
  mu_unif >= ci[1] && mu_unif <= ci[2]
})
prob_ci_contains_mu_unif <- mean(ci_contains_mu_unif)
cat("Uniform: Probability CI contains mu:", prob_ci_contains_mu_unif, "\n")

Uniform: Probability CI contains mu: 0.9347 

reject_H0_twosided_unif <- replicate(n_sim, {
  x <- runif(n_samp, min = 0, max = 1)
  t.test(x, mu = mu_unif, alternative = "two.sided")$p.value < 0.05
})
prob_reject_H0_twosided_unif <- mean(reject_H0_twosided_unif)
cat("Uniform: Probability of rejecting H0 (two-sided):", prob_reject_H0_twosided_unif, "\n")

Uniform: Probability of rejecting H0 (two-sided): 0.0672 

reject_H0_greater_unif <- replicate(n_sim, {
  x <- runif(n_samp, min = 0, max = 1)
  t.test(x, mu = mu_unif, alternative = "greater")$p.value < 0.05
})
prob_reject_H0_greater_unif <- mean(reject_H0_greater_unif)
cat("Uniform: Probability of rejecting H0 (greater):", prob_reject_H0_greater_unif, "\n")

Uniform: Probability of rejecting H0 (greater): 0.0535 

Conclusion: The distribution of \(X\) affects the accuracy of statistical inference procedures. When \(X\) is normal, coverage and error rates match theoretical expectations. When \(X\) is non-normal, these procedures may be less reliable, especially with moderate sample sizes.

3. Now let’s simulate with \(X\) following a normal distribution of \(N(1,4)\) (variance of 4 and standard devivation of 2).

(a) With a sample size of 10, do a hypothesis test for \(X\) with \(H_0: \mu = 0\) and \(H_a: \mu > 0\) with 10,000 simulations and \(\alpha = 0.05\). Find the \(\beta\) from your simulation. What is the estimated power of the test?

mu_alt <- 1
sd_alt <- 2
n_samp_alt <- 10

reject_H0_alt <- replicate(n_sim, {
  x <- rnorm(n_samp_alt, mean = mu_alt, sd = sd_alt)
  t.test(x, mu = 0, alternative = "greater")$p.value < 0.05
})
power_est <- mean(reject_H0_alt)
beta_est <- 1 - power_est
cat("Power (n=10, alpha=0.05):", power_est, "Beta:", beta_est, "\n")

Power (n=10, alpha=0.05): 0.4263 Beta: 0.5737 

(b) Change \(\alpha\) to 0.01 in (a) and find \(\beta\) and the power of test again.

reject_H0_alt_01 <- replicate(n_sim, {
  x <- rnorm(n_samp_alt, mean = mu_alt, sd = sd_alt)
  t.test(x, mu = 0, alternative = "greater")$p.value < 0.01
})
power_est_01 <- mean(reject_H0_alt_01)
beta_est_01 <- 1 - power_est_01
cat("Power (n=10, alpha=0.01):", power_est_01, "Beta:", beta_est_01, "\n")

Power (n=10, alpha=0.01): 0.1618 Beta: 0.8382 

(c) Change the sample size to 1000 and redo (a) and (b) again.

n_samp_alt_large <- 1000

reject_H0_alt_large <- replicate(n_sim, {
  x <- rnorm(n_samp_alt_large, mean = mu_alt, sd = sd_alt)
  t.test(x, mu = 0, alternative = "greater")$p.value < 0.05
})
power_est_large <- mean(reject_H0_alt_large)
beta_est_large <- 1 - power_est_large
cat("Power (n=1000, alpha=0.05):", power_est_large, "Beta:", beta_est_large, "\n")

Power (n=1000, alpha=0.05): 1 Beta: 0 

reject_H0_alt_large_01 <- replicate(n_sim, {
  x <- rnorm(n_samp_alt_large, mean = mu_alt, sd = sd_alt)
  t.test(x, mu = 0, alternative = "greater")$p.value < 0.01
})
power_est_large_01 <- mean(reject_H0_alt_large_01)
beta_est_large_01 <- 1 - power_est_large_01
cat("Power (n=1000, alpha=0.01):", power_est_large_01, "Beta:", beta_est_large_01, "\n")

Power (n=1000, alpha=0.01): 1 Beta: 0 

(d) How does this experiment verify the theory that we learn? Summarize your findings properly.

Power increases with larger samples and higher significance levels, while Type II error decreases under the same conditions. This demonstrates the importance of adequate sample size and appropriate choice of \(\alpha\) in hypothesis testing.