Monte Carlo Methods
Statistical Inference
Monte Carlo Simulations
Use Monte Carlo Simulations to show that transforming a data set to the T-Statistic will yield a t-distribution. Use \(H_0: \mu = 3\) and simulate 24 RV’s from \(N(3, 9)\). Testing how the distribution will look like assuming that the null hypothesis is true.
Answer:
Code
t_stat <- function(i){
x <- rnorm(24, 3, 3)
tt <- (mean(x) - 3)/(sd(x)/sqrt(24))
return(tt)
}
t_dist <- sapply(1:100000, t_stat)
xe <- seq(-5, 5, length.out = 100)
t_dist |> tibble(x = _) |>
ggplot(aes(x, y = ..density..)) +
geom_histogram() +
geom_line(data = tibble(x = xe, y = dt(xe, 23)),
mapping = aes(x,y),
col = "red", lwd = 2)
MC Sampling Distribution: \(\bar X\)
Code
sim1 <- function(i){
x <- rnorm(1000, mean = 5, sd = sqrt(2))
return(mean(x))
}
results <- replicate(10000, sim1(1))
xe <- seq(4.8, 5.2, length.out = 100)
tibble(x = results) |>
ggplot(aes(x, y = ..density..)) +
geom_histogram() +
geom_line(data = tibble(x = xe, y = dnorm(xe, 5, sqrt(2/1000))),
mapping = aes(x,y),
col = "red", lwd = 2)
MC Sampling Distribution: \(\hat s^2\)
Code
sim2 <- function(i){
x <- rnorm(1000, mean = 5, sd = sqrt(2))
return(var(x))
}
results <- replicate(10000, sim2(1))
xe <- seq(850, 1150, length.out = 100)
tibble(x = 999 * results / 2) |>
ggplot(aes(x, y = ..density..)) +
geom_histogram() +
geom_line(data = tibble(x = xe, y = dchisq(xe, 999)),
mapping = aes(x,y),
col = "red", lwd = 2)
Maximum
\[ f_{max} (x) = n f(x) F(x)^{n-1} \]
Code
sim3 <- function(i){
x <- rnorm(1000, mean = 5, sd = sqrt(2))
return(max(x))
}
results <- replicate(10000, sim3(1))
xe <- seq(8, 12, length.out = 100)
tibble(x = results) |>
ggplot(aes(x, y = ..density..)) +
geom_histogram() +
geom_line(data = tibble(x = xe,
y = 1000 * dnorm(xe, 5, sqrt(2)) * pnorm(xe, 5, sqrt(2))^999),
mapping = aes(x,y),
col = "red", lwd = 2)
Minimum
\[ f_{min}(x) = n\{1-F(x)\}^{n-1}f(x) \]
Code
sim4 <- function(i){
x <- rnorm(1000, mean = 5, sd = sqrt(2))
return(min(x))
}
results <- replicate(10000, sim4(1))
xe <- seq(-2.5, 2, length.out = 100)
tibble(x = results) |>
ggplot(aes(x, y = ..density..)) +
geom_histogram() +
geom_line(data = tibble(x = xe,
y = 1000 * (1 - pnorm(xe, 5, sqrt(2)))^999 * dnorm(xe, 5, sqrt(2))),
mapping = aes(x,y),
col = "red", lwd = 2)
t-test Example
Use a Monte Carlo Hypothesis Test for \(H_0: \mu = 3\) and \(H_a: \mu \neq 3\). Simulate fake data from 24 RV’s from \(N(3, 9)\). Testing how the distribution will look like assuming that the null hypothesis is true.
Answer:
Code
tstat <- (mean(x) - 3)/(sd(x)/sqrt(24))
t_stat_sim <- function(i){
x <- rnorm(24, 3, 3)
tt <- (mean(x) - 3)/(sd(x)/sqrt(24))
return(tt)
}
t_dist <- sapply(1:100000, t_stat_sim)
t_dist |> tibble(x = _) |>
ggplot(aes(x, y = ..density..)) +
geom_histogram() +
geom_vline(aes(xintercept = tstat), col = "red")#> [1] 69032#> [1] 0.6903231
Type 1 and 2 Error
Code
a1 <- data.frame(
x = c(2.5,3.5),
y = c(41, 41),
label = c("Yes", "No")
)
a2 <- data.frame(
x = c(1.9),
y = c(25, 35),
label = c("False", "True")
)
a3 <- tibble(x = c(1.7,3),
y = c(30, 43),
label = c("H0", "Reject H0"))
yay <- tibble(x = c(2.5, 3.5),
y = c(25, 35),
label = c("Yay!", "Yay!"))
type1 <- tibble(x = c(2.5, 3.5),
y = c(36.5, 26.5),
label = c("Type I Error", "Type II Error"))
type2 <- tibble(x = c(3.5, 2.5),
y = c(24.5, 34.5),
label = c("beta", "alpha"))
# basic graph
p <- ggplot() + theme_void()
# Add rectangles
p + annotate("rect",
xmin=c(2,3), xmax=c(3,4),
ymin=c(20,20), ymax=c(30,30),
alpha=0.2, color="green", fill="green") +
annotate("rect",
xmin=c(2,3), xmax=c(3,4),
ymin=c(30,30), ymax=c(40,40),
alpha=0.2, color="red", fill="red") +
annotate("rect",
xmin=c(3,2), xmax=c(4,3),
ymin=c(30, 20), ymax=c(40, 30),
alpha=0.8, color="royalblue", fill="royalblue") +
geom_text(data=a1, aes(x=x, y=y, label=label),
size=10 , fontface="bold" ) +
geom_text(data=a3, aes(x=x, y=y, label=label),
size=10 , fontface="bold.italic" ) +
geom_text(data=a2, aes(x=x, y=y, label=label),
size=10, angle = 90, fontface="bold" ) +
geom_text(data=yay, aes(x=x, y=y, label=label),
size=10, fontface="bold.italic" ) +
geom_text(data=type1, aes(x=x, y=y, label=label),
size=10, fontface="bold") +
geom_text(data=type2, aes(x=x, y=y, label=label),
size=10, fontface="bold", parse = T ) 
Type 1 and 2 Error
Code
x <- seq(-4, 4, length.out = 1000)
xx <- seq(-1, 7, length.out = 1000)
dt_one<-function(x){
y <- dnorm(x)
y[x < qnorm(.94)] <-NA
return(y)
}
dt_two<-function(x){
y <- dnorm(x, mean = 3)
y[x > 1.645] <-NA
return(y)
}
df1 <- tibble(x = x, y = dnorm(x))
df2 <- tibble(x = xx, y = dnorm(xx, mean = 3))
a1 <- tibble(x = c(0,3), y = c(0.43, 0.43), label = c("H0", "H1"))
a2 <- tibble(x = c(2.05,1.25), y = c(0.02, 0.02), label = c(paste("alpha"), "beta"))
df1 |>
ggplot(aes(x, y)) +
stat_function(fun = dt_one, geom = "area", fill = "red") +
geom_line() +
geom_line(aes(x, y), df2) +
stat_function(fun = dt_two, geom = "area", fill = "green") +
geom_vline(aes(xintercept = 1.645), lwd = 1.5) +
geom_text(data=a1, aes(x=x,y=y, label=label),
size=14 , fontface="bold.italic" ) +
geom_text(data=a2, aes(x=x,y=y, label=label),
size=10 , fontface="bold.italic", parse = T) +
theme_bw()
Conceptial
Code
x <- seq(-4, 4, length.out = 1000)
xx <- seq(-1, 7, length.out = 1000)
dt_one<-function(x){
y <- dnorm(x)
y[x < qnorm(.94)] <-NA
return(y)
}
dt_two<-function(x){
y <- dnorm(x, mean = 3)
y[x > 1.645] <-NA
return(y)
}
df1 <- tibble(x = x, y = dnorm(x))
df2 <- tibble(x = xx, y = dnorm(xx, mean = 3))
a1 <- tibble(x = c(0,3), y = c(0.43, 0.43), label = c("H0", "H1"))
a2 <- tibble(x = c(2.05,1.25), y = c(0.02, 0.02), label = c(paste("alpha"), "beta"))
df1 |>
ggplot(aes(x, y)) +
stat_function(fun = dt_one, geom = "area", fill = "red") +
geom_line() +
geom_line(aes(x, y), df2) +
stat_function(fun = dt_two, geom = "area", fill = "green") +
geom_vline(aes(xintercept = 1.645), lwd = 1.5) +
geom_text(data=a1, aes(x=x,y=y, label=label),
size=14 , fontface="bold.italic" ) +
geom_text(data=a2, aes(x=x,y=y, label=label),
size=10 , fontface="bold.italic", parse = T) +
theme_bw()