Running a Simulation Study

What is this?

• Already had a lecture on Monte Carlo Simulation Basics
• Recall the “General Recipe” for Monte Carlo:
  1. Write a function to carry out the experiment once.
  2. Use iteration to call that function a large number of times.
  3. Store and summarize results; set seed for replicability.
• Today: more elaborate examples that call for a more complicated recipe and more powerful tools.

🐱 More on purrr

• We used purr::map() in an earlier lecture. More today!
• References:
  • purrr cheatsheet
  • purrr tutorial
  • R for Data Science, Chapter 21
  • Advanced R Chapter 9
  • Functional Programming

What is a functional?

• A functional is a function that:
  • takes another function as an input or
  • returns another function as its output.
• Example: purrr::map(.x, .f) \[ \texttt{map}(\mathbf{x}, f) = \texttt{map}\left(\begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix},\, f\right) = \begin{bmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{bmatrix} \]

map(.x, .f) replaces for()

\[ \texttt{map}(\mathbf{x}, f) = \texttt{map}\left(\begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix},\, f\right) = \begin{bmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{bmatrix} \]

• .x is a list or vector; .f is a function
• map(x, f) gives same output as this for() loop:

results <- vector('list', length(x)) # pre-allocate empty list
for (i in seq_along(x)) { # safer than 1:length(x)
  results[[i]] <- f(x[[i]]) # x is a list; items could be *anything* 
} 

purrr:map() Example

sum_prod <- \(v) {
  # Return the sum and product of a vector v as a named vector
  c('sum' = sum(v), 'product' = prod(v))
}

x <- list(c(1, 1, 1), c(1, 2, 3), c(3, 3, 3))

library(tidyverse) # includes purrr
map(x, sum_prod)

[[1]]
    sum product 
      3       1 

[[2]]
    sum product 
      6       6 

[[3]]
    sum product 
      9      27 

Nice! But we may want to output something besides a list…

map_dfc() and map_dfr()

• map_dfc() binds columns (c), returns a data frame (df)

map_dfc(x, sum_prod) # looks for row names but there are none!

# A tibble: 2 × 3
   ...1  ...2  ...3
  <dbl> <dbl> <dbl>
1     3     6     9
2     1     6    27

• map_dfr() binds rows (r), returns a data frame (df)

map_dfr(x, sum_prod) # preserves column names

# A tibble: 3 × 2
    sum product
  <dbl>   <dbl>
1     3       1
2     6       6
3     9      27

Supplying a vector to map()

Here 1:5 is interpreted as though it were a list:

map(1:4, \(x) x^2)

[[1]]
[1] 1

[[2]]
[1] 4

[[3]]
[1] 9

[[4]]
[1] 16

Same works with map_dfc(), map_dfr() etc.

map_dbl(.x, .f) and friends

• Common special case:
  • .f returns a scalar
  • Ouput an atomic vector with same length as .x
• Specify the type of the output:
  • map_dbl() returns a double vector
  • map_chr() returns a character vector
  • map_int() returns an integer vector
  • map_lgl() returns a logical vector

map_xxx() Examples

map_dbl(1:10, \(x) x^2)

 [1]   1   4   9  16  25  36  49  64  81 100

map_int(c(5L, 10L, 13L, 7L, 3L), \(x) x %% 3)

[1] 2 1 1 1 0

map_chr(c('CORE', 'ERM', 'RULES'), tolower)

[1] "core"  "erm"   "rules"

map_lgl(c(TRUE, FALSE, TRUE, TRUE), \(x) !x)

[1] FALSE  TRUE FALSE FALSE

purrr::map2(.x, .y, .f)

\[ \texttt{map2}\left(\begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix},\, \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}, \, f\right) = \begin{bmatrix} f(x_1, y_1) \\ f(x_2, y_2) \\ \vdots \\ f(x_n, y_n) \end{bmatrix} \]

• .x and .y are vectors or lists of the same length
• .f is a function that takes two arguments
• Returns a list by default; variants map2_dfr(), map2_dfc(), map2_dbl() are analogous to map()

Parallel map: pmap(.l, .f)

\[ \texttt{pmap}\left(\begin{bmatrix} l_{11} & l_{12} & \cdots & l_{1r} \\ l_{21} & l_{22} & \cdots & l_{2r}\\ \vdots \\ l_{n1} & l_{n2} & \cdots & l_{nr} \end{bmatrix}, \, f\right) = \begin{bmatrix} f(l_{11}, l_{12}, \dots, l_{1r}) \\ f(l_{21}, l_{22}, \dots , l_{2r}) \\ \vdots \\ f(l_{n1}, l_{n2}, \dots, l_{nr}) \end{bmatrix} \]

• .l is a list of lists or a list of vectors (data frame)
• .f is a function that takes \(r\) arguments
• Returns a list by default; variants pmap_dfr(), pmap_dfc(), pmap_dbl() are analogous to map()

map2() and pmap() Examples

Argument names are arbitrary:

map2_dbl(6:10, 1:5, \(foo, bar) foo^bar)

[1]      6     49    512   6561 100000

df <- tibble(col1 = 1:5, col2 = 5:1, col3 = -2:2)
df

# A tibble: 5 × 3
   col1  col2  col3
  <int> <int> <int>
1     1     5    -2
2     2     4    -1
3     3     3     0
4     4     2     1
5     5     1     2

Argument names need to match df

pmap_dbl(df, \(col1, col2, col3) col3 / (col1 + col2)^2)

[1] -0.05555556 -0.02777778  0.00000000  0.02777778  0.05555556

💪 Exercise A - (3 min)

Recall this code snippet from our earlier lecture on simulation:

dice_sum <- \() {
# Roll a pair of fair, six-sided dice and return their sum
  die1 <- sample(1:6, 1)
  die2 <- sample(1:6, 1)
  die1 + die2
}
sims <- map_dbl(1:10000, \(i) dice_sum())

1. Use your new-found knowledge of purrr to explain line 7.
2. Write a for loop to replace line 7.

Sim Example - Estimators of \(\sigma^2\)

• \(X_1, ..., X_n \sim \text{iid N}(\mu, \sigma^2)\)
• \(S^2\) is the “usual” estimator; \(\hat{\sigma}^2\) is the MLE \[ S^2 \equiv \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X}_n)^2, \quad \widehat{\sigma}^2 \equiv \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X}_n)^2. \]
• \(S^2\) is unbiased; \(\mathbb{E}[\widehat{\sigma}^2] = (n-1)\sigma^2/n\).
• \(\text{Bias}(\widehat{\sigma}^2) = \mathbb{E}(\hat{\sigma}^2 - \sigma^2) = -\sigma^2/n\)

The Whole Kit & Caboodle

1. Write a function to generate simulation data.
2. Write a function to calculate your estimates.
3. Run the simulation for fixed parameters. Repeat many times:
   1. Draw simulation data.
   2. Calculate estimates.
4. Repeat step 4 over a grid range of parameter values.
5. Store and summarize the results.

1 - Function to Generate Sim Data

draw_sim_data <- \(n, s_sq) {
  rnorm(n, sd = sqrt(s_sq))
}

• For simplicity, take \(\mu = 0\) in our sim (doesn’t matter).
• Allow n and s_sq to vary

2 - Function to Calculate Estimates

get_estimates <- \(x) {
  c('usual' = var(x),
    'MLE' = mean((x - mean(x))^2)) # divide by n; not (n - 1)
}

possibly()

• Sometimes an estimator may fail to exist for some datasets.
• If an error occurs, don’t want the simulation to crash.
• possibly() a function to fail gracefully, suppressing errors and storing a designated value in case of failure:

x <- list('Core ERM', 3, 4)
map_dbl(x, log)

Error in `map_dbl()`:
ℹ In index: 1.
Caused by error:
! non-numeric argument to mathematical function

map_dbl(x, possibly(log, NA_real_))

[1]       NA 1.098612 1.386294

3 - Run Sim for Fixed Parameters

# Generate list of three simulation datasets as an example
set.seed(1693)
sim_datasets <- map(1:3, \(i) draw_sim_data(n = 5, s_sq = 1))
sim_datasets

[[1]]
[1] -0.5070343 -0.8750463 -0.1883818 -1.1772701  1.8960400

[[2]]
[1]  0.3824529  1.0291447 -1.3998332  1.3983002 -0.3950949

[[3]]
[1]  0.1204153 -0.9611377 -0.3870477  0.7757843  0.6520975

# Compute the estimates for all three datasets; bind rows
map_dfr(sim_datasets, get_estimates)

# A tibble: 3 × 2
  usual   MLE
  <dbl> <dbl>
1 1.47  1.18 
2 1.27  1.01 
3 0.527 0.421

3 - Run Sim for Fixed Parameters

Now that we have the idea, do more replications!

set.seed(1693)
nreps <- 5000
sim_datasets <- map(1:nreps, \(i) draw_sim_data(n = 5, s_sq = 1))
sim_estimates <- map_dfr(sim_datasets, get_estimates)
sim_estimates

# A tibble: 5,000 × 2
   usual   MLE
   <dbl> <dbl>
 1 1.47  1.18 
 2 1.27  1.01 
 3 0.527 0.421
 4 0.858 0.686
 5 0.863 0.690
 6 0.942 0.754
 7 0.521 0.417
 8 2.84  2.27 
 9 0.941 0.753
10 0.775 0.620
# ℹ 4,990 more rows

3 - Run Sim for Fixed Parameters

Using sim_estimates we can approximate the bias of \(\widehat{\sigma}^2\)

# Sim parameters: n = 5,  s_sq = 1, 
n <- 5
s_sq <- 1
sim_estimates |> 
  summarize(bias_true = -s_sq / n, bias_sim = mean(MLE - s_sq))

# A tibble: 1 × 2
  bias_true bias_sim
      <dbl>    <dbl>
1      -0.2   -0.197

Great: our simulation appears to work!

4 - Run Sim over Parameter Grid

Wrap step 3 from above into a function:

run_sim <- \(n, s_sq, nreps = 5000) { # default: nreps = 5000
  map(1:nreps, \(i) draw_sim_data(n, s_sq)) |> 
    map_dfr(get_estimates)
}
set.seed(1693)
run_sim(n = 5, s_sq = 1)

# A tibble: 5,000 × 2
   usual   MLE
   <dbl> <dbl>
 1 1.47  1.18 
 2 1.27  1.01 
 3 0.527 0.421
 4 0.858 0.686
 5 0.863 0.690
 6 0.942 0.754
 7 0.521 0.417
 8 2.84  2.27 
 9 0.941 0.753
10 0.775 0.620
# ℹ 4,990 more rows

4 - Run Sim over Parameter Grid

Use expand_grid() from tidyr to set up parameter grid

sim_params <- expand_grid(n = 3:5, s_sq = 1:3)
sim_params

# A tibble: 9 × 2
      n  s_sq
  <int> <int>
1     3     1
2     3     2
3     3     3
4     4     1
5     4     2
6     4     3
7     5     1
8     5     2
9     5     3

5 - Summarize Results

get_summary_stats <- function(df) {
  c('usual_mean' = mean(df$usual), 
    'MLE_mean' = mean(df$MLE), 
    'usual_var' = var(df$usual), 
    'MLE_var' = var(df$MLE))
}
summary_stats <- map_dfr(sim_results, get_summary_stats)
summary_stats <- bind_cols(sim_params, summary_stats)
summary_stats

# A tibble: 9 × 6
      n  s_sq usual_mean MLE_mean usual_var MLE_var
  <int> <int>      <dbl>    <dbl>     <dbl>   <dbl>
1     3     1      0.990    0.660     0.966   0.429
2     3     2      1.96     1.31      3.82    1.70 
3     3     3      2.93     1.95      8.58    3.81 
4     4     1      0.995    0.746     0.688   0.387
5     4     2      2.00     1.50      2.60    1.46 
6     4     3      2.95     2.21      5.80    3.26 
7     5     1      0.997    0.797     0.484   0.310
8     5     2      2.02     1.62      2.09    1.34 
9     5     3      3.03     2.42      4.71    3.01 

5 - Summarize Results

summary_stats |> 
  mutate(MLE_bias = -s_sq / n, 
         MLE_sim_bias = MLE_mean - s_sq,
         usual_sim_bias = usual_mean - s_sq,
         MLE_rmse = sqrt(MLE_sim_bias^2 + MLE_var),
         usual_rmse = sqrt(usual_sim_bias^2 + usual_var)) |> 
  select(MLE_bias, MLE_sim_bias, MLE_rmse, usual_rmse)

# A tibble: 9 × 4
  MLE_bias MLE_sim_bias MLE_rmse usual_rmse
     <dbl>        <dbl>    <dbl>      <dbl>
1   -0.333       -0.340    0.738      0.983
2   -0.667       -0.691    1.48       1.96 
3   -1           -1.05     2.22       2.93 
4   -0.25        -0.254    0.672      0.829
5   -0.5         -0.498    1.31       1.61 
6   -0.75        -0.789    1.97       2.41 
7   -0.2         -0.203    0.592      0.696
8   -0.4         -0.383    1.22       1.45 
9   -0.6         -0.578    1.83       2.17 

All Roads Lead to Rome

• Many different ways to set up the same simulation study.
• Trade-offs: speed / memory usage / ease of testing
• I could have written a function to run the simulation once:

draw_sim_rep <- \(...) {
  draw_sim_data() |> 
    get_estimates()
}
get_estimates2 <- \(nreps = 5000, ...) {
  map_dbl(1:nreps, \(i) draw_sim_rep(...))
}

• Could also have written the code for draw_sim_data() and get_estimates() within draw_sim_data()

Kahneman - Thinking Fast & Slow

Amos [Tversky] and his students Tom Gilovich and Robert Vallone caused a stir with their study of misperceptions of randomness in basketball. The “fact” that players occasionally acquire a hot hand is generally accepted by players, coaches, and fans. The inference is irresistible: a player sinks three or four baskets in a row and you cannot help forming the causal judgment that this player is now hot, with a temporarily increased propensity to score. Players on both teams adapt to this judgment—teammates are more likely to pass to the hot scorer and the defense is more likely to doubleteam.

Kahneman - Thinking Fast & Slow

Analysis of thousands of sequences of shots led to a disappointing conclusion: there is no such thing as a hot hand in professional basketball, either in shooting from the field or scoring from the foul line. Of course, some players are more accurate than others, but the sequence of successes and missed shots satisfies all tests of randomness. The hot hand is entirely in the eye of the beholders, who are consistently too quick to perceive order and causality in randomness. The hot hand is a massive and widespread cognitive illusion.

Oops!

Simulation is more powerful than derivation. Don’t be clever; be right.