Simulation & Variability — Summarize simulated results and interpret variability (intro) - R Programming

SECTION 1

Historical Context & Motivation

Long before digital computers existed, scientists and mathematicians grappled with problems that had no closed-form analytical solution — integrals over irregular domains, stochastic system behaviors, and combinatorial explosions that rendered enumeration impractical. The breakthrough insight was deceptively simple: if you can simulate a process many times and record the outcomes, you can approximate any quantity of interest — its mean, its spread, and its full distribution — without ever solving the equations directly. This idea underpins the entire family of techniques now known as Monte Carlo methods, and understanding the variability in simulated results is the key to extracting trustworthy conclusions from them.

1777 • Buffon's Needle Problem

Comte de Buffon proposed estimating π by dropping a needle on parallel lines — one of the earliest probabilistic simulation arguments, laying conceptual groundwork for Monte Carlo reasoning.

1946 • Monte Carlo Method at Los Alamos

Stanislaw Ulam and John von Neumann formalized Monte Carlo simulation for neutron diffusion calculations during the Manhattan Project, exploiting the ENIAC computer's ability to run thousands of random trials.

1976 • Efron's Bootstrap

Bradley Efron introduced the bootstrap, a resampling-based simulation technique that uses the empirical distribution of observed data to estimate sampling variability without parametric assumptions.

1993 • R Language Created

Ross Ihaka and Robert Gentleman developed R at the University of Auckland, providing an open-source environment tailor-made for statistical computing, including vectorized random-number generation and rich plotting facilities.

2000s–Present • Simulation-Based Inference in CS & Data Science

Monte Carlo simulations are now routine in machine learning (cross-validation, hyperparameter tuning), systems research (queuing models, network simulations), and Bayesian computation (MCMC), making simulation literacy essential for computer scientists.

The central question this lesson addresses is straightforward yet profound: once you have generated a large collection of simulated outcomes in R, how do you summarize those results, and how do you interpret the variability you observe? Answering this requires understanding summary statistics, distributional shapes, the role of the number of replications, and the connection between simulation variability and inferential uncertainty.

SECTION 2

Core Principles & Definitions

Before diving into R code, it is important to establish a precise vocabulary. Simulation in the statistical computing sense means using a pseudo-random number generator (PRNG) to produce synthetic data according to a specified probabilistic model, then studying the behavior of some statistic or quantity computed from that data. Each independent run of this process is called a replication, and the collection of values obtained across all replications constitutes the simulated sampling distribution. Variability in this distribution is not noise to be ignored — it is the fundamental quantity that allows us to make probabilistic statements about the real world.

Replication

A single execution of the simulation procedure. In R, one iteration inside a loop or one call to replicate() that produces a scalar or vector result.

Summary Statistics

Numerical summaries that condense the simulated distribution: mean(), sd(), quantiles, and range. These estimate the center, spread, and shape of the underlying process.

Simulation Variability

The spread in results across replications caused by randomness. Larger samples per replication reduce within-replication noise; more replications reduce Monte Carlo error in the summary.

Monte Carlo Error

The uncertainty in a simulation-based estimate attributable to using a finite number of replications. It shrinks proportionally to 1/√n where n is the replication count.

Reproducibility via set.seed()

R's set.seed() fixes the PRNG state, ensuring identical results across runs — critical for debugging, peer review, and scientific reproducibility.

✦

KEY TAKEAWAY

Think of a simulation study like running the same experiment in thousands of parallel universes. Each universe gives a slightly different answer because of random chance. The spread of answers across universes is the variability, the average answer is your best estimate, and running more universes makes that average more precise — like reducing the margin of error in a poll by surveying more people.

SECTION 3

Visual Explanation — Anatomy of a Simulation

The following diagram illustrates the complete simulation workflow in R: you define a data-generating process, run it many times via replicate(), collect the results into a vector, and then compute summary statistics and visualize the distribution. The key insight is that each replication yields a different value due to randomness, and the histogram of all those values reveals both the central tendency and the variability of the simulated quantity.

The top row shows the four-step workflow: define a data-generating process, replicate it many times, summarize the results, and interpret. The histogram below shows a typical simulated sampling distribution of the sample mean from 10,000 replications. The dashed pink line marks the grand mean, and the amber bracket indicates ±1 standard deviation — the simulation variability.

Notice that the histogram is approximately normal, which is expected by the Central Limit Theorem when each replication computes a sample mean. The spread of this histogram — quantified by its standard deviation — is the simulation-estimated standard error. If the population standard deviation is σ = 15 and each sample has n = 30 observations, the theoretical standard error is σ/√n = 15/√30 ≈ 2.74, matching the simulation result closely. This agreement is the hallmark of a well-designed simulation study.

SECTION 4

Mathematical Framework

The mathematical backbone of simulation-based analysis rests on a few key results from probability theory. Understanding these equations lets you predict how precise a simulation estimate will be before you even run it, and lets you verify your code by comparing simulated values to theoretical ones.

MONTE CARLO ESTIMATOR

θ̂ = (1/B) × Σᵢ₌₁ᴮ T(Xᵢ*)

Where θ̂ is the simulation estimate of the parameter θ, B is the number of replications, Xᵢ* is the simulated dataset in replication i, and T(·) is the statistic of interest (e.g., mean()).

STANDARD ERROR OF THE MEAN

SE(X̄) = σ / √n

The theoretical standard error of the sample mean where σ is the population standard deviation and n is the sample size within each replication. This governs within-replication variability.

MONTE CARLO STANDARD ERROR

MCSE(θ̂) = s_T / √B

Where $s_T$ is the sample standard deviation of the B computed statistics and B is the number of replications. This quantifies the Monte Carlo error — the imprecision due to running a finite simulation. Quadrupling B halves the MCSE.

SIMULATION-BASED CONFIDENCE INTERVAL

CI = [θ̂ − z* × MCSE, θ̂ + z* × MCSE]

For a 95% CI, z* ≈ 1.96. Alternatively, use the 2.5th and 97.5th percentiles of the simulated distribution directly via quantile(results, c(0.025, 0.975)).

⚠

Two Sources of Variability

It is critical to distinguish sampling variability (inherent randomness in each sample, governed by SE = σ/√n) from Monte Carlo variability (imprecision from using finite replications, governed by MCSE =

s_T

/√B). Increasing the sample size n within each replication reduces the first; increasing the number of replications B reduces the second.

SECTION 5

Detailed Breakdown — R Functions for Simulation

R's vectorized architecture and built-in random number generators make it exceptionally well-suited for Monte Carlo work. The key functions form a coherent pipeline: generate random data with distribution functions like rnorm(), rbinom(), or sample(); automate repetition with replicate(); and summarize with standard statistical functions. The diagram below maps R functions to each stage of the workflow.

Three panels show the R functions for each stage: random generation (left), repetition (center), and summarization (right). The bottom panel shows a complete working example that generates 10,000 sample means.

R Function	Purpose	Returns
`set.seed(s)`	Fix PRNG state for reproducibility	NULL (side effect)
`rnorm(n, mean, sd)`	Generate n normal random variates	Numeric vector of length n
`sample(x, n, replace)`	Draw n items from vector x	Vector of length n
`replicate(B, expr)`	Evaluate expr B times independently	Vector (or matrix) of length B
`quantile(x, probs)`	Compute specified percentiles	Named numeric vector
`hist(x, breaks)`	Plot histogram of simulated values	Histogram plot (invisible list)

Key R functions for simulation workflows

SECTION 6

Worked Example — Estimating a Proportion via Simulation

Suppose a CS professor claims that 35% of students prefer Python over R. You want to simulate the sampling distribution of the sample proportion p̂ when surveying n = 50 students, using B = 5,000 replications. You will then summarize the results and interpret the variability.

Simulating the Sampling Distribution of p̂

1
Step 1 — Set up the simulation parameters

Define the true proportion p = 0.35, sample size n = 50, and number of replications B = 5000. Set the random seed for reproducibility: set.seed(123); p <- 0.35; n <- 50; B <- 5000

Parameters: p = 0.35, n = 50, B = 5,000
2
Step 2 — Define the single-replication function

Each replication draws n = 50 Bernoulli outcomes from rbinom(1, size = 50, prob = 0.35) and divides by n to obtain p̂. Equivalently, we can write: sim_prop <- function() { sum(rbinom(50, 1, 0.35)) / 50 }. Note that rbinom(1, 50, 0.35) / 50 is a more vectorized alternative.

Each call returns a single p̂ value between 0 and 1
3
Step 3 — Run all replications

Use p_hats <- replicate(B, rbinom(1, n, p) / n) to produce a numeric vector of 5,000 simulated sample proportions. This single line replaces a for-loop and is idiomatic R.

p_hats is a numeric vector of length 5,000
4
Step 4 — Compute summary statistics

Calculate the mean, standard deviation, and key quantiles: mean(p_hats) ≈ 0.3503, sd(p_hats) ≈ 0.0676, and quantile(p_hats, c(0.025, 0.975)) ≈ [0.22, 0.48]. The theoretical SE is √(p(1−p)/n) = √(0.35 × 0.65 / 50) ≈ 0.0675, confirming the simulation.

Simulated mean ≈ 0.3503, SD ≈ 0.0676, 95% interval ≈ [0.22, 0.48]
5
Step 5 — Interpret the variability

The standard deviation of the 5,000 p̂ values (≈ 0.068) estimates the standard error of the sample proportion. This tells us that in repeated surveys of 50 students, the observed proportion would typically differ from the true 0.35 by about 6.8 percentage points. The 95% simulation interval [0.22, 0.48] indicates that it would be unsurprising to observe p̂ values anywhere in that range — important context if a single survey yields p̂ = 0.28 and you wonder whether the true p might still be 0.35.

Variability is substantial: surveys of 50 students cannot pin down a proportion to better than about ±7 percentage points

SECTION 7

Strengths & Limitations of Simulation-Based Analysis

Simulation is extraordinarily flexible — it can tackle problems for which no closed-form solution exists, and it provides tangible, visual evidence of variability. However, like any computational technique, it has trade-offs that practitioners must understand.

Strengths	Limitations
Works for any statistic — no need to derive formulas for SE analytically	Results are approximate: always subject to Monte Carlo error
Provides the full distribution, not just point estimates	Computationally expensive for very large B or complex models
Makes abstract concepts (CLT, SE, sampling distributions) concrete and visual	Garbage in, garbage out: results depend on correctly specified models
Easily extended to multi-step or correlated data processes	PRNG quality matters; rare-event simulation may need specialized methods
Reproducible via set.seed() — crucial for collaborative CS projects	Different seeds yield different results; must report Monte Carlo uncertainty

Strengths vs. limitations of Monte Carlo simulation in R

✦

KEY TAKEAWAY

Simulation is to statistics what unit testing is to software engineering: it does not prove correctness, but it gives you strong empirical evidence about how a system behaves under realistic conditions. Just as you would not ship software with a single test case, you should not run a simulation with too few replications — more runs reduce Monte Carlo error and increase confidence in your conclusions.

SECTION 8

Connection to Advanced Theory

The introductory simulation techniques covered in this lesson serve as the foundation for a rich family of advanced methods. As you progress through statistical computing, you will encounter scenarios where the basic replicate() pattern extends into significantly more powerful frameworks. Understanding how variability behaves in simple simulations prepares you to reason about convergence, efficiency, and reliability in these more complex settings.

This Lesson (Intro)	Advanced Extension	Key New Concept
Simulate from known distributions (rnorm, rbinom)	Bootstrap: resample from observed data	Empirical distribution replaces parametric assumption
Fixed B replications	Adaptive stopping rules	Monitor MCSE and stop when precision target is met
Independent replications	Markov Chain Monte Carlo (MCMC)	Correlated samples; convergence diagnostics required
Single statistic per replication	Permutation tests	Simulate null distribution to compute p-values
Summary via mean/sd	Power analysis via simulation	Proportion of replications rejecting H₀ estimates power

From intro simulation to advanced statistical computing

The thread connecting all of these advanced techniques back to this lesson is the principle that variability across replications is informative, not a nuisance. In the bootstrap, the variability of resampled statistics estimates the sampling distribution. In MCMC, the variability of the chain after burn-in characterizes the posterior distribution. In permutation tests, the variability of the test statistic under random permutations defines the null distribution. Mastering the basics of summarizing and interpreting simulation variability here gives you a transferable mental model for all of these methods.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain the difference between sampling variability (the standard error of a statistic) and Monte Carlo variability (Monte Carlo standard error). If you run a simulation with B = 10,000 replications, each drawing n = 100 observations from N(50, 10²), which source of variability does sd(sim_means) estimate, and which does sd(sim_means) / sqrt(B) estimate?

Show answer

The value sd(sim_means) estimates the standard error of the sample mean, which is σ/√n = 10/√100 = 1.0. This is the sampling variability — the inherent spread you would see if you repeated the study infinitely many times. The quantity sd(sim_means) / sqrt(B) ≈ 1.0/√10000 = 0.01 estimates the Monte Carlo standard error: how precisely our simulation-based estimate of that SE itself is pinned down. In short, the first measures the real-world phenomenon; the second measures the imprecision introduced by the finite simulation.

Question 2: BASIC CALCULATION

You simulate B = 4,000 replications, each computing the median of n = 20 observations from Uniform(0, 1). You find that mean(sim_medians) = 0.4988 and sd(sim_medians) = 0.0878. (a) What is the Monte Carlo standard error of your estimate of the true mean of the sampling distribution? (b) Construct an approximate 95% confidence interval for the true expected median.

Show answer

(a) The MCSE is sd(sim_medians)/√B = 0.0878/√4000 ≈ 0.0878/63.25 ≈ 0.00139. (b) A 95% CI is 0.4988 ± 1.96 × 0.00139 = [0.4961, 0.5015]. Since the theoretical median of U(0,1) is exactly 0.5 and this interval contains 0.5, the simulation is consistent with theory. This illustrates how even B = 4,000 replications yield a very precise estimate of the expected value, while the SD of 0.088 captures the considerable replication-to-replication variability in the sample median.

Question 3: INTERMEDIATE

Write R code to estimate the probability that the maximum of 5 independent Uniform(0,1) random variables exceeds 0.95. Use B = 10,000 replications. Then compute the simulation-based estimate and a 95% CI for the true probability. Hint: the theoretical answer is 1 − 0.95⁵.

Show answer

Code:

set.seed(7); B <- 10000; results <- replicate(B, max(runif(5)) > 0.95); p_hat <- mean(results); se <- sqrt(p_hat * (1 - p_hat) / B)

. The theoretical probability is 1 − 0.95⁵ = 1 − 0.7738 = 0.2262. The simulation should yield p̂ ≈ 0.226 with SE ≈ √(0.226 × 0.774 / 10000) ≈ 0.00418, giving a 95% CI of approximately [0.218, 0.234]. Here, results is a logical vector; mean(results) coerces TRUEs to 1s and computes the proportion. The SE formula uses the Bernoulli variance formula p(1−p)/B since each replication is a binary outcome.

Question 4: APPLIED

A load balancer distributes incoming requests to 3 servers. Each server's processing time is Exponential with rate λ = 2 (mean = 0.5 seconds). A request is assigned to a server uniformly at random, and the total time is the processing time of the chosen server. Suppose 10 requests arrive simultaneously and are distributed independently. Write an R simulation with B = 5,000 replications to estimate the mean and standard deviation of the total time to process all 10 requests (i.e., the maximum processing time across all 10). Summarize and interpret the variability.

Show answer

Code:

set.seed(99); B <- 5000; max_times <- replicate(B, { times <- rexp(10, rate = 2); max(times) }); cat('Mean:', mean(max_times), '\nSD:', sd(max_times), '\n95% interval:', quantile(max_times, c(0.025, 0.975)))

. With 10 i.i.d. Exp(2) variables, the expected maximum is H₁₀/2 where H₁₀ is the 10th harmonic number ≈ 2.928/2 ≈ 1.464. The simulation should yield a mean around 1.46 seconds and SD around 0.50 seconds. The 95% interval [≈ 0.67, ≈ 2.58] shows substantial variability — a system designer cannot assume worst-case processing will consistently be near the mean. This variability is critical for capacity planning in distributed systems.

Question 5: CRITICAL THINKING

A colleague runs a simulation with B = 100 replications and reports that the mean of their simulated statistic is 42.3 with SD = 8.1. They conclude that 'the true value is about 42.' (a) Compute the Monte Carlo standard error and the 95% CI for the true mean. (b) Now suppose they increase to B = 10,000. By what factor does the MCSE shrink? (c) Critically evaluate whether B = 100 is sufficient. Under what circumstances might B = 100 be acceptable, and when would it be dangerously misleading?

Show answer

(a) MCSE = 8.1/√100 = 0.81, so the 95% CI is 42.3 ± 1.96 × 0.81 = [40.71, 43.89]. This is a wide interval — the true expected value could plausibly be anywhere from about 41 to 44. (b) At B = 10,000, MCSE = 8.1/√10000 = 0.081, which is a factor of 10 smaller (since √(10000/100) = 10). The 95% CI narrows to [42.14, 42.46]. (c) B = 100 might be acceptable for a rough exploratory check or when each replication is very expensive (e.g., training a deep learning model). However, it is dangerously misleading when precise estimates are needed, when the underlying distribution is heavy-tailed (the SD of 8.1 might itself be poorly estimated with only 100 values), or when the simulation result will inform a high-stakes decision. A good practice is to always report the MCSE alongside any simulation-based estimate, analogous to reporting confidence intervals with any statistical estimate.

SUMMARY

Lesson Summary

Summary

This lesson introduced the foundational workflow for Monte Carlo simulation in R: define a data-generating process using functions like rnorm() or rbinom(), automate repetition with replicate(), summarize results using mean(), sd(), and quantile(), and visualize the simulated sampling distribution with hist(). The standard deviation across replications estimates the standard error of the statistic, while the Monte Carlo standard error (MCSE = s/√B) quantifies the imprecision from using a finite number of replications.

The critical takeaway is that variability in simulated results is the signal, not noise — it reveals the inherent uncertainty in the process being studied. By distinguishing between sampling variability (controlled by sample size n) and Monte Carlo variability (controlled by replication count B), you gain the ability to design simulations that are both informative and computationally efficient. Always use set.seed() for reproducibility and always report the MCSE alongside your simulation estimates. These foundational skills transfer directly to advanced techniques including the bootstrap, permutation tests, and MCMC.