Lab 2 Solutions
# function to flip a coin with probability p
flip <- function(p) {
sample(c(0,1), size=1, prob=(c(1-p,p)))
}
# test out flip function
flip(0.5)
#> [1] 0
# function to generate a sample of size n
generate_sample <- function(n,p) {
Y <- c()
for (i in 1:n) {
Y[i] <- flip(p)
}
Y
}
# test out generate_sample function
generate_sample(10,0.5)
#> [1] 1 0 0 1 1 1 1 1 0 0
# carry out monte carlo simulations
n <- 10
p <- 0.5
nsims <- 1000 # need to pick large number of monte carlo simulations
mc_est <- c() # vector to hold estimation results
mc_var <- c() # vector to hold estimated variance
for (i in 1:nsims) {
Y <- generate_sample(n,p)
mc_est[i] <- mean(Y)
mc_var[i] <- var(Y)
}
# compute bias
bias <- mean(mc_est) - p
bias
#> [1] 0.0046
# compute sampling variance
var <- var(mc_est)
var
#> [1] 0.0240629
# compute mean squared error
mse <- bias^2 + var
mse
#> [1] 0.02408406
H0 <- p
t <- sqrt(n)*(mc_est - H0) / sqrt(mc_var)
ggplot(data.frame(t=t), aes(x=t)) +
geom_histogram(bins=30) +
theme_bw()
#> Warning: Removed 2 rows containing non-finite values
#> (`stat_bin()`).
rej <- mean(1*(abs(t) >= 1.96))
rej
#> [1] 0.107
# since we are going to do this over and over, let's write a function to do it
mc_sim <- function(n, p, H0) {
mc_est <- c() # vector to hold estimation results
mc_var <- c() # vector to hold estimated variance
for (i in 1:nsims) {
Y <- generate_sample(n,p)
mc_est[i] <- mean(Y)
mc_var[i] <- var(Y)
}
# compute bias
bias <- mean(mc_est) - p
# compute sampling variance
var <- var(mc_est)
# compute mean squared error
mse <- bias^2 + var
t <- sqrt(n)*(mc_est - H0) / sqrt(mc_var)
hist_plot <- ggplot(data.frame(t=t), aes(x=t)) +
geom_histogram(bins=30) +
theme_bw()
rej <- mean(1*(abs(t) >= 1.96))
# print results
print(paste0("bias: ", round(bias,4)))
print(paste0("var : ", round(var,4)))
print(paste0("mse : ", round(mse,4)))
print(hist_plot)
print(paste0("rej : ", round(rej,4)))
}
mc_sim(50, 0.5, 0.5)
#> [1] "bias: 0.0014"
#> [1] "var : 0.0056"
#> [1] "mse : 0.0056"
#> [1] "rej : 0.085"
mc_sim(50, 0.5, 0.6)
#> [1] "bias: -0.0015"
#> [1] "var : 0.0048"
#> [1] "mse : 0.0048"
#> [1] "rej : 0.339"
mc_sim(50, 0.5, 0.9)
#> [1] "bias: 0.0011"
#> [1] "var : 0.005"
#> [1] "mse : 0.005"
#> [1] "rej : 1"
mc_sim(1000, 0.5, 0.6)
#> [1] "bias: 4e-04"
#> [1] "var : 2e-04"
#> [1] "mse : 2e-04"
#> [1] "rej : 1"
mc_sim(10, 0.95, 0.95)
#> [1] "bias: -0.0058"
#> [1] "var : 0.005"
#> [1] "mse : 0.005"
#> Warning: Removed 555 rows containing non-finite values
#> (`stat_bin()`).
#> [1] "rej : 0.556"
mc_sim(50, 0.95, 0.95)
#> [1] "bias: -0.001"
#> [1] "var : 9e-04"
#> [1] "mse : 9e-04"
#> Warning: Removed 59 rows containing non-finite values
#> (`stat_bin()`).
#> [1] "rej : 0.064"
mc_sim(1000, 0.95, 0.95)
#> [1] "bias: 1e-04"
#> [1] "var : 0"
#> [1] "mse : 0"
#> [1] "rej : 0.057"