Quick summary
Summarize this blog with AI
Introduction
A Monte Carlo estimate approximates a defined quantity by repeated simulated trials. More trials reduce simulation error under the model; they do not validate the model’s assumptions. Here the target is the probability that two independent fair six-sided dice sum to at least 10. We can enumerate the exact answer, then use it to check the simulation.
Tested with R 4.5.1. Run the code blocks in order. Install the named packages before starting; inputs are supplied below.
Define the event and an exact reference
There are 36 equally likely ordered pairs and six successful pairs. This benchmark makes a programming mistake visible.
pairs <- expand.grid(first = 1:6, second = 1:6)
exact <- mean(pairs$first + pairs$second >= 10)
stopifnot(exact == 1/6)
RNGkind("Mersenne-Twister", "Inversion", "Rejection")
set.seed(2026)
B <- 20000L
first <- sample.int(6L, B, replace = TRUE)
second <- sample.int(6L, B, replace = TRUE)
totals <- first + second
hit <- totals >= 10
estimate <- mean(hit)
print(c(estimate = estimate, exact = exact))
stopifnot(length(hit) == B, all(totals >= 2 & totals <= 12))
Each simulated trial uses two fresh independent draws with replacement. The seed makes this particular run reproducible; it does not make the estimate exact. The expected estimate is 0.16395, compared with the exact probability 1/6, approximately 0.1667.
Quantify simulation error in the probability estimate
The indicator is Bernoulli. Its estimated Monte Carlo standard error is sqrt(p_hat×(1−p_hat)/B). An exact binomial interval applies to the simulated success count under these independent trials.
mcse <- sqrt(estimate * (1 - estimate) / B)
interval <- binom.test(sum(hit), B, conf.level = 0.95)$conf.int
print(round(c(estimate = estimate, mcse = mcse,
lower = interval[1], upper = interval[2]), 5))
stopifnot(mcse > 0, interval[1] < estimate, interval[2] > estimate)
stopifnot(abs(estimate - exact) < 5 * sqrt(exact * (1 - exact) / B))
The expected rounded output is estimate 0.16395, standard error 0.00262, and interval [0.15884, 0.16915]. The interval describes Monte Carlo uncertainty about the probability under the specified dice model. It is not uncertainty about whether physical dice are fair. The five-standard-error assertion is a regression check for this fixed fixture, not a claim that every random run must pass it.
Do not confuse an outcome range with an estimator interval
Quantiles of totals summarize individual dice outcomes. They answer a different question from the binomial interval above.
outcome_range <- quantile(totals, c(.025, .975), type = 1)
print(outcome_range) # 2 and 12 in this run
stopifnot(identical(unname(outcome_range), c(2, 12)))
checkpoints <- c(1000L, 5000L, 20000L)
progress <- data.frame(trials = checkpoints,
estimate = vapply(checkpoints, function(n) mean(hit[seq_len(n)]), numeric(1)),
theoretical_mcse = sqrt(exact * (1 - exact) / checkpoints))
print(progress)
stopifnot(all(diff(progress$theoretical_mcse) < 0))
Theoretical standard error falls as 1/sqrt(B). The observed estimate need not move closer to the truth at every checkpoint. To halve standard error, use approximately four times as many independent trials; this does not remove bias from a wrong simulation model.
Validate trial counts and handle rare events
Reject an empty simulation rather than reporting NaN as a probability. Zero observed successes does not prove a probability is zero.
check_trials <- function(n) {
if (!is.numeric(n) || length(n) != 1L || !is.finite(n) ||
n < 1 || n != floor(n) || n > 1e6) stop("Use an integer trial count from 1 to 1000000")
invisible(TRUE)
}
check_trials(B)
for (bad in list(0, numeric(0), NA_real_, Inf, 2.5, "100"))
stopifnot(inherits(tryCatch(check_trials(bad), error = identity), "error"))
zero_success_interval <- binom.test(0, 1000)$conf.int
stopifnot(zero_success_interval[1] == 0, zero_success_interval[2] > 0)
The one-million limit is a teaching budget, not an R limit. For larger studies, aggregate batches without storing every draw. Parallel random streams and importance sampling need their own validated designs; unrelated random draws are not automatically MCMC or importance sampling.
