4.18 Lab 3 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"