Chapter 7 Simulation and Reproducibility

Simulation helps test intuition, evaluate estimators, and understand sampling variation. A reproducible simulation specifies the random seed, data-generating process, number of repetitions, and quantity being estimated.

7.1 Distribution functions

R uses consistent prefixes for probability distributions:

  • d* returns a density or probability mass.
  • p* returns a cumulative probability.
  • q* returns a quantile.
  • r* generates random values.
dnorm(0)
#> [1] 0.3989423
pnorm(1.96)
#> [1] 0.9750021
qnorm(0.975)
#> [1] 1.959964

For a continuous distribution, dnorm(0) is a density height, not the probability that a continuous random variable equals exactly zero.

7.2 Reproducible random generation

set.seed() resets R’s random-number stream so the same code produces the same sequence.

set.seed(2026)
first_draw <- rnorm(5)

set.seed(2026)
second_draw <- rnorm(5)

identical(first_draw, second_draw)
#> [1] TRUE

Set a seed at the simulation boundary. Repeatedly resetting it inside a helper function can accidentally generate identical samples across calls.

7.3 Sampling

Sampling without replacement selects distinct values. Sampling with replacement permits a value to appear multiple times and underlies the nonparametric bootstrap.

set.seed(2026)
sample(1:10, size = 4, replace = FALSE)
#> [1] 9 1 6 5
sample(1:6, size = 8, replace = TRUE)
#> [1] 3 4 4 5 4 2 2 3

7.4 Simulating a linear model

Make the data-generating process explicit before fitting the model.

set.seed(2026)
simulation_size <- 500
predictor <- rnorm(simulation_size)
noise <- rnorm(simulation_size, mean = 0, sd = 2)
outcome <- 0.5 + 2 * predictor + noise

simulated_fit <- lm(outcome ~ predictor)
round(coef(simulated_fit), 2)
#> (Intercept)   predictor 
#>        0.48        1.96

The fitted coefficients will not equal the generating parameters exactly because each sample contains random variation.

7.5 Monte Carlo estimation

A Monte Carlo estimate is an average over repeated random experiments. Its simulation error generally decreases as the number of repetitions increases.

set.seed(2026)
repetitions <- 10000
dice_totals <- replicate(
  repetitions,
  sum(sample(1:6, size = 2, replace = TRUE))
)

mean(dice_totals >= 10)
#> [1] 0.1617

Report the seed and repetition count so another analyst can reproduce the result and assess simulation precision.

7.6 Bootstrap uncertainty

The nonparametric bootstrap repeatedly samples observed values with replacement and recomputes a statistic.

observed_revenue <- c(20, 25, 30, 35, 50, 70)

set.seed(2026)
bootstrap_means <- replicate(
  2000,
  mean(sample(observed_revenue, replace = TRUE))
)

round(quantile(bootstrap_means, c(0.025, 0.975)), 2)
#>  2.5% 97.5% 
#> 25.83 52.50

This interval reflects the empirical sample and bootstrap assumptions. It is not a substitute for checking whether observations are independent and representative of the target population.

7.7 Timing and profiling

Optimize only after identifying a real bottleneck. system.time() measures an expression, while Rprof() and dedicated profilers show where a program spends time.

system.time({
  result <- expensive_function(data)
})

Rprof("profile.out")
result <- expensive_function(data)
Rprof(NULL)
summaryRprof("profile.out")

Profiling examples are not executed here because they require a real workload and create output files.

7.8 Interview checkpoints

  • Explain the d, p, q, and r distribution prefixes.
  • Explain what set.seed() guarantees and what it does not.
  • Distinguish sampling with and without replacement.
  • Describe a data-generating process for a simple linear model.
  • Explain what uncertainty a nonparametric bootstrap represents.