Monte Carlo Simulation Basics in R - Solutions

Exercise A - (4 min)

  1. What happens if I don’t supply size to sample()?
  2. Does it matter if I set replace = TRUE or FALSE here?
sample(1:6, size = 1)
  1. Write a function called rbern() that uses sample() to make a sequence of \(n\) iid Bernoulli\((p)\) draws. Its arguments are n, the number of draws, and p, the probability of success. It should return a vector n of zeros and ones. Check your results by making a large number of draws and verifying that the fraction of ones is approximately \(p\).

Solutions

  1. It sets size to length(x)
  2. No: since we’re only making a single draw from 1:6, it doesn’t matter if we draw with or without replacement.
  3. E.g.
rbern <- \(n, p) {
  sample(0:1, size = n, replace = TRUE, prob = c(1 - p, p))
}
mean(rbern(1e6, 0.2)) - 0.2
[1] -0.000149
mean(rbern(1e6, 0.4)) - 0.4
[1] 0.00114
mean(rbern(1e6, 0.6)) - 0.6
[1] -9.1e-05
mean(rbern(1e6, 0.8)) - 0.8
[1] 0.000133

Exercise B - (10 min)

  1. Consult ?dexp() and then carry out the following:
    1. Plot the density and CDF of an Exponential(1) RV.
    2. Plot a histogram of 500 draws from an Exponential(1) distribution.
  2. Consult ?dpois() and then carry out the following:
    1. Plot the mass function and CDF of a Poisson(5) RV.
    2. Plot a barchart of 500 draws from a Poisson(5) distribution.
  3. Suppose that we observe a sample of \(n=16\) iid normal RVs with unknown mean \(\mu\) and known variance \(\sigma^2 = 1\).
    1. A confidence interval for \(\mu\) takes the form \(\bar{X}\pm \text{ME}\) where \(\text{ME}\) is the margin of error. Use qnorm() to compute \(\text{ME}\) for an 89% interval.
    2. Simulate 16 draws from a \(N(0,\sigma=4)\) distribution and construct the interval from part (a). Does it cover the true value of \(\mu\)?

Solution

library(tidyverse)

# Part 1a
exponential1 <- tibble(x = seq(0, 5, 0.01)) |> 
  mutate(density = dexp(x, 1), cdf = pexp(x, 1))

exponential1 |> 
  ggplot(aes(x, density)) + 
  geom_line() +
  labs(title = 'Exponential(1) Density')

exponential1 |> 
  ggplot(aes(x, cdf)) + 
  geom_line() +
  labs(title = 'Exponential(1) CDF')

# Part 1b
set.seed(690872)
tibble(x = rexp(500, 1)) |> 
  ggplot(aes(x)) +
  geom_histogram(bins = 25)

# Part 2a
# pmf
tibble(x = 0:20) |>  
  mutate(p = dpois(x, 5)) |> 
  ggplot(aes(x, p)) + 
  geom_point() +
  geom_segment(aes(xend = x, yend = 0)) +
  labs(title = 'Poisson(5) pmf', y =