sample(1:6, size = 1)
Monte Carlo Simulation Basics in R - Solutions
Exercise A - (4 min)
- What happens if I don’t supply
size
tosample()
? - Does it matter if I set
replace = TRUE
orFALSE
here?
- Write a function called
rbern()
that usessample()
to make a sequence of \(n\) iid Bernoulli\((p)\) draws. Its arguments aren
, the number of draws, andp
, the probability of success. It should return a vectorn
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
- It sets
size
tolength(x)
- No: since we’re only making a single draw from
1:6
, it doesn’t matter if we draw with or without replacement. - E.g.
<- \(n, p) {
rbern 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)
- Consult
?dexp()
and then carry out the following:- Plot the density and CDF of an Exponential(1) RV.
- Plot a histogram of 500 draws from an Exponential(1) distribution.
- Consult
?dpois()
and then carry out the following:- Plot the mass function and CDF of a Poisson(5) RV.
- Plot a barchart of 500 draws from a Poisson(5) distribution.
- Suppose that we observe a sample of \(n=16\) iid normal RVs with unknown mean \(\mu\) and known variance \(\sigma^2 = 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. - 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\)?
- A confidence interval for \(\mu\) takes the form \(\bar{X}\pm \text{ME}\) where \(\text{ME}\) is the margin of error. Use
Solution
library(tidyverse)
# Part 1a
<- tibble(x = seq(0, 5, 0.01)) |>
exponential1 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 =