knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'tidyr' was built under R version 4.4.3
## Warning: package 'purrr' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(stats)
set.seed(123)

Problem 1: Normal Distribution 𝑁(0,1)

Let’s say a random variable 𝑋 follows a normal distribution of 𝑁(0,1).

a. Draw 5 samples of 𝑋 and use 𝑥¯ 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 𝜇=0.

n_sims <- 10000
n <- 5
µ_true <- 0
sigma <- 1

cover_1a <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma/sqrt(n)
  lower <- xbar-qnorm(0.975)*se
  upper <- xbar+qnorm(0.975)*se
  (lower <= µ_true) && (upper >= µ_true)
})

p_cover_1a <- mean(cover_1a)
round(p_cover_1a, 4)
## [1] 0.9506

b. Do a hypothesis test for 𝑋 with 𝐻0:𝜇=0 and 𝐻𝑎:𝜇≠0 with 10,000 simulations and 𝛼=0.05. Find the probability of rejecting 𝐻0 by comparing the p-value and 𝛼.

n_sims <- 10000
n <- 5
µ_true <- 0
sigma <- 1
alpha <- 0.05

reject_two_1b <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 2 * (1 - pnorm(abs(z)))
  pval < alpha
})

reject_prob_two_1b <- mean(reject_two_1b)
round(reject_prob_two_1b, 4)
## [1] 0.0486

c. Redo (b) with 𝐻𝑎:𝜇>0.

n_sims <- 10000
n <- 5
µ_true <- 0
sigma <- 1
alpha <- 0.05

reject_right_1c <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

reject_prob_right_1c <- mean(reject_right_1c)
round(reject_prob_right_1c, 4)
## [1] 0.0503

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

n_sims <- 10000
n <- 100
µ_true <- 0
sigma <- 1
alpha <- 0.05

cover_1d <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma/sqrt(n)
  lower <- xbar-qnorm(0.975)*se
  upper <- xbar+qnorm(0.975)*se
  (lower <= µ_true) && (upper >= µ_true)
})

p_cover_1d <- mean(cover_1d)
round(p_cover_1d, 4)
## [1] 0.9479
n_sims <- 10000
n <- 100
µ_true <- 0
sigma <- 1
alpha <- 0.05

reject_two_1d <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 2 * (1 - pnorm(abs(z)))
  pval < alpha
})

reject_prob_two_1d <- mean(reject_two_1d)
round(reject_prob_two_1d, 4)
## [1] 0.0505
n_sims <- 10000
n <- 100
µ_true <- 0
sigma <- 1
alpha <- 0.05

reject_right_1d <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

reject_prob_right_1d <- mean(reject_right_1d)
round(reject_prob_right_1d, 4)
## [1] 0.0506

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

Answers: 1. The probability of the estimated interval is very close to 95%, the same as the confident level we choose 2. The probability of rejecting 𝐻0 is around 0.05, closed to the alpha we picked. 3. The sample size has limited influence on the results we already have.

Problem 2: Uniform Distribution Unif(-1,1)

Repeat all steps in Problem 1 with X following a uniform distribution Unif(-1,1). How does the distribution of X affect the results?

Unif(-1, 1) -> The ture µ is still 0, while sigma change to √(1/3)

2a.

n <- 5
alpha <- 0.05
µ_true <- 0
sigma <- sqrt(1/3)

cover <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma/sqrt(n)
  lower <- xbar-qnorm(0.975)*se
  upper <- xbar+qnorm(0.975)*se
  (lower <= µ_true) && (upper >= µ_true)
})

p_cover <- mean(cover)
round(p_cover, 4)
## [1] 0.9525

2b.

n <- 5
alpha <- 0.05
µ_true <- 0
sigma <- sqrt(1/3)

reject_two_2b <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 2 * (1 - pnorm(abs(z)))
  pval < alpha
})

reject_prob_two <- mean(reject_two_2b)
round(reject_prob_two, 4)
## [1] 0.0497

2c.

n <- 5
alpha <- 0.05
µ_true <- 0
sigma <- sqrt(1/3)

reject_right_2c <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

reject_prob_right <- mean(reject_right_2c)
round(reject_prob_right, 4)
## [1] 0.05

2d.

n <- 100
alpha <- 0.05
µ_true <- 0
sigma <- sqrt(1/3)

cover <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma/sqrt(n)
  lower <- xbar-qnorm(0.975)*se
  upper <- xbar+qnorm(0.975)*se
  (lower <= µ_true) && (upper >= µ_true)
})

p_cover <- mean(cover)
round(p_cover, 4)
## [1] 0.9519
n <- 100
alpha <- 0.05
µ_true <- 0
sigma <- sqrt(1/3)

reject_two_2d <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 2 * (1 - pnorm(abs(z)))
  pval < alpha
})

reject_prob_two <- mean(reject_two_2d)
round(reject_prob_two, 4)
## [1] 0.0553
n <- 100
alpha <- 0.05
µ_true <- 0
sigma <- sqrt(1/3)

reject_right_2d <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

reject_prob_right <- mean(reject_right_2d)
round(reject_prob_right, 4)
## [1] 0.0481

2e.

Answers: 1. The probability of the estimated interval is very close to 95%, the same as the confident level we choose 2. The probability of rejecting 𝐻0 is around 0.05, closed to the alpha we picked. 3. The sample size has limited influence on the results we already have.

Comments: The distribution of X doesn’t really affect the results.

Problem 3: Normal Distribution N(1,4)

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

µ_true <- 1
sigma <- 2

a. With a sample size of 10, do a hypothesis test for 𝑋 with 𝐻0:𝜇=0 and 𝐻𝑎:𝜇>0 with 10,000 simulations and 𝛼=0.05. Find the 𝛽 from your simulation. What is the estimated power of the test?

n <- 10
alpha <- 0.05

rejected_3a <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

power_3a <- mean(rejected_3a)
beta_3a <- 1 - power_3a
c(beta = round(beta_3a,4), power = round(power_3a,4))
##   beta  power 
## 0.5256 0.4744

b. Change 𝛼 to 0.01 in (a) and find 𝛽 and the power of test again.

n <- 10
alpha <- 0.01

rejected_3b <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

power_3b <- mean(rejected_3b)
beta_3b <- 1 - power_3b
c(beta = round(beta_3b,4), power = round(power_3b,4))
##  beta power 
## 0.772 0.228

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

n <- 1000
alpha <- 0.05

rejected_3ca <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

power_3ca <- mean(rejected_3ca)
beta_3ca <- 1 - power_3ca
c(beta = round(beta_3ca,4), power = round(power_3ca,4))
##  beta power 
##     0     1
n <- 1000
alpha <- 0.01

rejected_3cb <- replicate(n_sims, {
  x <- rnorm(n, mean = µ_true, sd = sigma)
  xbar <- mean(x)
  se <- sigma / sqrt(n)
  z <- (xbar - 0) / se
  pval <- 1 - pnorm(z)
  pval < alpha
})

power_3cb <- mean(rejected_3cb)
beta_3cb <- 1 - power_3cb
c(beta = round(beta_3cb,4), power = round(power_3cb,4))
##  beta power 
##     0     1

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

Answers: 1. When we choose alpha at 0.05, the probability of having Type II error (Fail to reject wrong Null Hypothesis) is around 50-50 2. As the alpha getting lower (from 0.05 to 0.01), the probability of having Type II error has significantly increase. (Beta increase) 3. As the sample size increase, the probability of having Type II error will decrease.