In the previous chapter, we used computation to estimate parameters from data. Here we go in the other direction: given a probability model, how do we generate observations from it?
We start with uniform random numbers, use them to construct one-dimensional distributions, and then build mixtures and random vectors. At each step, we ask the same three questions: What is the construction? Why does it have the desired distribution? How can we check it?
Generating random variables is the first step in a simulation study. Once we can generate data with a known truth, we can study an estimator’s bias, a confidence interval’s coverage, or a test’s power.
The main methods in this chapter follow a common pattern:
| Available building block | Construction | Target |
|---|---|---|
| Uniform probability levels | Apply a quantile function | A specified univariate distribution |
| An easy proposal distribution | Accept or reject proposals | A distribution with a known density |
| Familiar independent variables | Transform or combine them | Related distributions, sums, and mixtures |
| Independent standard normals | Apply a matrix or normalize lengths | Correlated vectors and random directions |
A pseudorandom number generator produces a deterministic sequence from an initial state. A well-designed generator makes that sequence useful as an approximation to independent random draws. In the mathematical arguments below, we work with ideal independent uniforms; software approximates them with finite precision.
In R, use set.seed() to reproduce an example. The same seed, generator settings, software environment, and sequence of calls reproduce the same simulation. The seed identifies a computation; it does not certify the probability model.
set.seed(8670)
u_first <- runif(5)
set.seed(8670)
u_repeated <- runif(5)
rbind(first = u_first, repeated = u_repeated)
#> [,1] [,2] [,3] [,4] [,5]
#> first 0.3114911 0.7752026 0.874258 0.7551533 0.9447904
#> repeated 0.3114911 0.7752026 0.874258 0.7551533 0.9447904
identical(u_first, u_repeated)
#> [1] TRUESet the seed before a simulation experiment, then let the generator advance. Resetting the same seed inside every replication creates repeated copies of the same random input.
For a reproducible analysis, also record RNGkind() and sessionInfo(). For parallel work, use a framework that manages separate random-number streams, such as the streams supported by R’s parallel package. See the R documentation on random-number generation.
Sampling with replacement allows an observation to appear more than once. Sampling without replacement makes later draws depend on which observations have already been selected.
set.seed(8670)
observed_waits <- c(2, 3, 3, 5, 7, 12) # illustrative waiting times
sample(observed_waits, size = 6, replace = TRUE)
#> [1] 12 5 3 12 12 3
sample(observed_waits, size = 6, replace = FALSE)
#> [1] 12 3 2 7 3 5The empirical distribution places probability \(1/n\) on each of the \(n\) observed records, counting repeated values separately. Drawing \(n\) times independently from this distribution is the basic nonparametric bootstrap resampling step. Sampling all \(n\) records without replacement only permutes them, so the sample mean cannot change.
This is our first example of an important principle: the sampling rule is part of the statistical model.
R uses four related prefixes. For example, dnorm() evaluates a normal density, pnorm() evaluates its CDF, qnorm() evaluates its quantile function, and rnorm() generates draws.
| Distribution | Density or PMF | CDF | Quantile | Generator | Parameters to check |
|---|---|---|---|---|---|
| Uniform | dunif |
punif |
qunif |
runif |
min, max |
| Normal | dnorm |
pnorm |
qnorm |
rnorm |
mean, sd |
| Exponential | dexp |
pexp |
qexp |
rexp |
rate |
| Gamma | dgamma |
pgamma |
qgamma |
rgamma |
shape, rate or scale |
| Beta | dbeta |
pbeta |
qbeta |
rbeta |
shape1, shape2 |
| Binomial | dbinom |
pbinom |
qbinom |
rbinom |
size, prob |
| Geometric | dgeom |
pgeom |
qgeom |
rgeom |
prob; counts failures |
| Poisson | dpois |
ppois |
qpois |
rpois |
lambda |
| Negative binomial | dnbinom |
pnbinom |
qnbinom |
rnbinom |
size, prob or mu |
| Chi-square | dchisq |
pchisq |
qchisq |
rchisq |
df |
| Student’s \(t\) | dt |
pt |
qt |
rt |
df |
| \(F\) | df |
pf |
qf |
rf |
df1, df2 |
| Lognormal | dlnorm |
plnorm |
qlnorm |
rlnorm |
meanlog, sdlog |
See R’s distribution reference for the available families.
rnorm(n, sd = 4) has variance \(16\). An exponential rate of \(4\) gives mean \(1/4\). In the gamma family, scale = 1 / rate. A lognormal meanlog is the mean of \(\log X\), not the mean of \(X\).
A generator can run without errors while representing the wrong model if its parameterization is misunderstood.
Imagine dividing the interval \((0,1)\) into probability levels. A uniform draw selects a level without favoring one part of that interval. The quantile function converts that probability level into a value on the scale of the target distribution.
For any CDF \(F\), define
\[ F^{-1}(u)=\inf\{x\in\mathbb{R}:F(x)\geq u\}, \qquad 0<u<1. \]
It is the first value of \(x\) at which accumulated probability reaches the level \(u\).
If \(U\sim\operatorname{Unif}(0,1)\), then \(X=F^{-1}(U)\) has CDF \(F\). Independent uniforms transformed separately give independent observations from \(F\).
Proof. The generalized inverse satisfies
\[ \begin{aligned} P\{F^{-1}(U)\leq x\} &=P\{U\leq F(x)\}\\ &=F(x). \end{aligned} \]
Thus the generated values have the required CDF.
The condition is \(F(x)\geq u\), not \(F(x)=u\). A discrete CDF jumps over some probability levels. The generalized inverse still assigns every \(0<u<1\) to an outcome.
If \(X\) has a continuous CDF \(F\), then \(F(X)\sim\operatorname{Unif}(0,1)\). This is the reverse direction of the inverse-transform construction.
Under a correctly specified continuous predictive distribution, transformed observations have a uniform marginal distribution. Dependence must be checked separately. If \(F\) is estimated from the same observations, formal tests need to account for that estimation.
For discrete \(X\), \(F(X)\) is not uniform. An optional randomized version spreads each probability mass over its corresponding interval:
\[ U=F(X^-)+V\{F(X)-F(X^-)\}, \]
where \(V\) is an independent uniform and \(F(X^-)\) is the CDF just below \(X\).
Input: a target CDF \(F\) and sample size \(n\). Output: \(n\) independent draws from \(F\).
Consider
\[ f(x)=3x^2,\qquad 0<x<1. \]
First integrate the density, then solve for the observation:
\[ F(x)=\int_0^x 3t^2\,dt=x^3, \qquad u=x^3\ \Longrightarrow\ x=u^{1/3}. \]
Thus \(X=U^{1/3}\). For one draw, \(u=0.125\) gives \(x=0.5\); checking \(F(0.5)=0.125\) confirms the direction of the transformation.
u <- 0.125 # a given probability level
x <- u^(1 / 3) # convert it to an observation
c(x = x, F_x = x^3)
#> x F_x
#> 0.500 0.125Why do the observations concentrate near 1? Equal-width intervals near 1 contain more probability because the CDF rises faster there. The inverse maps more uniform probability levels into those intervals.
The arrows in Figure 4.2 show what the inverse does: choose a height \(u\), move across to the CDF, then read the corresponding \(x\) on the horizontal axis.
quantile_paths <- data.frame(
u = c(0.125, 0.5, 0.875),
level = factor(c("u = 0.125", "u = 0.500", "u = 0.875"))
)
quantile_paths$x <- quantile_paths$u^(1 / 3)
quantile_grid <- data.frame(x = seq(0, 1, length.out = 401))
ggplot(quantile_grid, aes(x, x^3)) +
geom_line(color = "grey25", linewidth = 1) +
geom_segment(
data = quantile_paths, inherit.aes = FALSE,
aes(x = 0, xend = x, y = u, yend = u, color = level),
linewidth = 0.8, linetype = "dashed", show.legend = FALSE,
arrow = grid::arrow(length = grid::unit(0.12, "inches"))
) +
geom_segment(
data = quantile_paths, inherit.aes = FALSE,
aes(x = x, xend = x, y = u, yend = 0, color = level),
linewidth = 0.8, linetype = "dashed", show.legend = FALSE,
arrow = grid::arrow(length = grid::unit(0.12, "inches"))
) +
geom_point(
data = quantile_paths, inherit.aes = FALSE,
aes(x, u, color = level), size = 2.5
) +
annotate("text", x = 0.28, y = 0.97, label = "F(x) = x^3", size = 4.5) +
scale_color_manual(values = c("#0072B2", "#009E73", "#D55E00")) +
scale_x_continuous(breaks = seq(0, 1, 0.25)) +
scale_y_continuous(breaks = seq(0, 1, 0.25)) +
coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
labs(x = "Generated value x", y = "Cumulative probability F(x)", color = NULL) +
theme(legend.position = "bottom")
set.seed(8670)
n <- 5000
u <- runif(n)
x_power <- u^(1 / 3)
plot_data <- rbind(
data.frame(value = u, distribution = "1. Uniform input"),
data.frame(value = x_power, distribution = "2. Transformed output")
)
grid <- seq(0, 1, length.out = 301)
curve_data <- rbind(
data.frame(value = grid, density = 1,
distribution = "1. Uniform input"),
data.frame(value = grid, density = 3 * grid^2,
distribution = "2. Transformed output")
)
ggplot(plot_data, aes(value)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
boundary = 0, fill = "lightblue", color = "white") +
geom_line(data = curve_data, aes(y = density),
color = "firebrick", linewidth = 1) +
facet_wrap(~ distribution) +
labs(x = "Value", y = "Density")
Here \(\mathbb{E}(X)=3/4\) and \(\operatorname{Var}(X)=3/80\). A simulation provides an independent numerical check.
knitr::kable(booktabs = TRUE,
data.frame(
quantity = c("Mean", "Variance"),
simulated = c(mean(x_power), var(x_power)),
theoretical = c(3 / 4, 3 / 80)
),
digits = 4
)| quantity | simulated | theoretical |
|---|---|---|
| Mean | 0.7531 | 0.7500 |
| Variance | 0.0372 | 0.0375 |
Agreement should be approximate. Random samples are not required to reproduce population moments exactly.
For \(T\sim\operatorname{Exp}(\lambda)\) with rate \(\lambda>0\),
\[ F(t)=1-e^{-\lambda t}, \qquad F^{-1}(u)=-\frac{\log(1-u)}{\lambda}. \]
Hence we can generate waiting times using
\[ T=-\frac{\log(1-U)}{\lambda}. \]
Since \(U\) and \(1-U\) have the same distribution, \(-\log(U)/\lambda\) also works. Both formulas have the same marginal distribution.
set.seed(8670)
lambda <- 0.7
u <- runif(10000)
waiting_time <- -log1p(-u) / lambda
c(
simulated_mean = mean(waiting_time),
theoretical_mean = 1 / lambda,
simulated_prob_above_3 = mean(waiting_time > 3),
theoretical_prob_above_3 = exp(-lambda * 3)
)
#> simulated_mean theoretical_mean simulated_prob_above_3
#> 1.4497408 1.4285714 0.1247000
#> theoretical_prob_above_3
#> 0.1224564log1p(-u) computes \(\log(1-u)\) accurately when \(u\) is close to zero. This connects random-number generation to the numerical stability issues discussed earlier.
The formulas \(-\log(U)/\lambda\) and \(-\log(1-U)/\lambda\) have the same marginal distribution. Applying both to the same \(U\) produces dependent outputs. Use independent uniforms when independence is required.
Exponential waiting times model the time between events in a homogeneous Poisson process. The constant rate is a substantive assumption: it implies a memoryless waiting time. Equipment aging or time-varying arrival rates may require a different distribution.
More generally, if a nonnegative event time has continuous cumulative hazard \(H(t)\) and survival function \(S(t)=e^{-H(t)}\) tending to zero, then
\[ T=H^{-1}\{-\log(U)\}. \]
This turns survival-model simulation into an inversion problem. For example, the Weibull cumulative hazard \(H(t)=(t/b)^a\) gives \(T=b\{-\log(U)\}^{1/a}\).
Exercise 4.1 (In Class: Rate or Mean? (2 minutes)) Multiple choice. To generate exponential waiting times with mean 5, use (A) rexp(n, rate = 5) or (B) rexp(n, rate = 1 / 5).
If the rate doubles, does the mean waiting time double or halve? Explain in one sentence.
Solution 4.1.
rate = 1 / 5, so the correct choice is B.An inverse may be inconvenient algebraically even when the CDF is easy to evaluate. For each \(u\), solve
\[ F(x)-u=0. \]
For example, consider
\[ F(x)=\frac{x+x^3}{2},\quad 0\leq x\leq 1, \qquad f(x)=\frac{1+3x^2}{2}. \]
The CDF is strictly increasing, so each \(u\in(0,1)\) has exactly one root in \((0,1)\). We can use the bracketing methods from the optimization chapter.
custom_cdf <- function(x) (x + x^3) / 2 # used only on [0, 1]
custom_quantile <- function(u) {
vapply(u, function(ui) {
uniroot(
function(x) custom_cdf(x) - ui,
interval = c(0, 1), tol = 1e-10
)$root
}, numeric(1))
}
set.seed(8670)
u <- runif(2000)
x_numeric <- custom_quantile(u)
c(
largest_cdf_residual = max(abs(custom_cdf(x_numeric) - u)),
simulated_mean = mean(x_numeric),
theoretical_mean = 5 / 8
)
#> largest_cdf_residual simulated_mean theoretical_mean
#> 4.630307e-11 6.305761e-01 6.250000e-01There are two distinct errors: the root solver introduces numerical error, while finite simulation introduces Monte Carlo error. Tightening tol reduces the former; increasing the number of independent draws reduces the latter. Neither repairs a wrong CDF.
Suppose only values in \((a,b)\) are eligible for a study. If \(F\) is continuous and \(F(b)>F(a)\), then the conditional CDF within this interval is
\[ P(X\leq x\mid a<X<b) =\frac{F(x)-F(a)}{F(b)-F(a)}. \]
Consequently,
\[ X=F^{-1}\left[F(a)+U\{F(b)-F(a)\}\right] \]
generates from the truncated distribution.
set.seed(8670)
a <- 1
b <- 2
u <- runif(5000)
truncated_normal <- qnorm(pnorm(a) + u * (pnorm(b) - pnorm(a)))
c(
simulated_mean = mean(truncated_normal),
theoretical_mean = (dnorm(a) - dnorm(b)) / (pnorm(b) - pnorm(a)),
fraction_in_interval = mean(truncated_normal > a & truncated_normal < b)
)
#> simulated_mean theoretical_mean fraction_in_interval
#> 1.388024 1.383169 1.000000Truncation changes which individuals enter the dataset. Censoring instead records incomplete information about some individuals who remain in the dataset, such as knowing that a failure time exceeds the study duration. These require different simulation rules. The direct CDF subtraction above is suitable for this moderate interval; extreme tails require numerically stable tail probabilities or a specialized truncated-distribution generator.
For ordered support values \(x_1<x_2<\cdots\), the quantile rule selects the first value whose cumulative probability reaches \(u\):
\[ X=x_j \quad\text{when}\quad F(x_{j-1})<U\leq F(x_j). \]
Each outcome receives a subinterval of \((0,1)\) whose length equals its probability.
For \(X\sim\operatorname{Bernoulli}(p)\), the generalized inverse gives
\[ X=\mathbf{1}\{U>1-p\}. \]
The threshold divides the unit interval into pieces of lengths \(1-p\) and \(p\).
The same idea applies to three categories with probabilities \(0.2\), \(0.5\), and \(0.3\). In Figure 4.4, interval length represents probability. A uniform draw at \(u=0.62\) selects the middle category because \(0.2<0.62\leq0.7\).
probability_strip <- data.frame(
left = c(0, 0.2, 0.7), right = c(0.2, 0.7, 1),
category = factor(c("Low", "Medium", "High"),
levels = c("Low", "Medium", "High")),
label = c("Low\np = 0.2", "Medium\np = 0.5", "High\np = 0.3")
)
ggplot(probability_strip) +
geom_rect(aes(xmin = left, xmax = right, ymin = 0.15, ymax = 0.7,
fill = category), color = "white", linewidth = 1.5) +
geom_text(aes(x = (left + right) / 2, y = 0.425, label = label),
size = 4.5, lineheight = 1.2) +
annotate("segment", x = 0.62, xend = 0.62, y = 1.02, yend = 0.73,
linewidth = 0.9,
arrow = grid::arrow(length = grid::unit(0.12, "inches"))) +
annotate("text", x = 0.62, y = 1.15, label = "U = 0.62 selects Medium",
size = 4.3) +
scale_fill_manual(values = c("#A9D4E5", "#A8DAB5", "#F4C49B"),
guide = "none") +
scale_x_continuous(breaks = c(0, 0.2, 0.7, 1), expand = expansion(mult = 0.04)) +
coord_cartesian(xlim = c(0, 1), ylim = c(0, 1.28)) +
labs(x = "Uniform probability level", y = NULL) +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
panel.grid = element_blank())
set.seed(8670)
p <- 0.4
bernoulli <- as.integer(runif(10000) > 1 - p)
labels <- c("Low", "Medium", "High")
probabilities <- c(0.2, 0.5, 0.3)
cutoffs <- cumsum(probabilities)
u <- runif(10000)
# First cutoff that reaches u: this also handles u at a cutoff.
index <- vapply(u, function(ui) which(ui <= cutoffs)[1], integer(1))
category <- factor(labels[index], levels = labels)
mean(bernoulli)
#> [1] 0.4094
prop.table(table(category))
#> category
#> Low Medium High
#> 0.1966 0.5008 0.3026
# Built-in alternative for the same categorical distribution.
head(sample(labels, size = 10, replace = TRUE, prob = probabilities))
#> [1] "Medium" "High" "High" "Medium" "Low" "High"For the geometric distribution, use R’s convention: \(X\) counts failures before the first success. Then
\[ P(X=x)=p(1-p)^x,\qquad x=0,1,2,\ldots, \]
and \(F(x)=1-(1-p)^{x+1}\). Solving \(F(x-1)<u\leq F(x)\) gives
\[ F^{-1}(u)= \left\lceil\frac{\log(1-u)}{\log(1-p)}\right\rceil-1. \]
set.seed(8670)
p <- 0.25
u <- runif(10000)
failures <- ceiling(log1p(-u) / log1p(-p)) - 1
c(
simulated_mean = mean(failures),
theoretical_mean = (1 - p) / p,
largest_difference_from_qgeom = max(abs(failures - qgeom(u, prob = p)))
)
#> simulated_mean theoretical_mean
#> 3.0503 3.0000
#> largest_difference_from_qgeom
#> 0.0000The number of trials including the success is \(X+1\), with mean \(1/p\). Mixing up these conventions creates an off-by-one modeling error. See R’s geometric distribution documentation.
For a Poisson distribution, cumulative probabilities can be constructed from
\[ p_0=e^{-\lambda},\qquad p_{k+1}=\frac{\lambda}{k+1}p_k, \qquad F(k+1)=F(k)+p_{k+1}. \]
Search for the first \(k\) satisfying \(F(k)\geq u\). This explains an inversion algorithm, but the recurrence can underflow for large \(\lambda\). In practice, use qpois(runif(n), lambda) or the specialized generator rpois(n, lambda). An r function need not implement the simple algorithm used in our derivation.
The inverse-transform method is especially convenient when a quantile is available. When a density is easier to work with than its CDF, the next method gives another route.
Inverse transformation requires a usable quantile function. Acceptance-rejection requires a different resource: an easy proposal distribution that can cover the target.
Let \(f\) be the normalized target density and \(g\) a proposal density. We need a finite constant \(M\) such that
\[ f(y)\leq M g(y) \]
everywhere on the target support. In particular, \(g(y)\) must be positive wherever \(f(y)>0\).
Input: target density \(f\), proposal density \(g\), a valid bound \(M\), and sample size \(n\).
Let \(f\) and \(g\) be normalized densities with \(f(y)\leq M g(y)\) throughout the target support. In the algorithm above:
Proof. For continuous proposals,
\[ P(\text{accept}) =\int g(y)\frac{f(y)}{M g(y)}\,dy =\frac{1}{M}. \]
The density of an accepted proposal is therefore
\[ f_{Y\mid\text{accept}}(y) =\frac{g(y)P(\text{accept}\mid Y=y)} {P(\text{accept})} =f(y). \]
The proposal visits some regions too frequently. The acceptance rule thins those visits by a location-specific amount. After thinning, the relative frequencies match the target.
For discrete distributions, use sums instead of integrals in the proof. The number of proposals through one acceptance is geometric on \(\{1,2,\ldots\}\) with mean \(M\). The realized cost varies from run to run.
This produces independent accepted draws, unlike a Markov chain method whose successive states are generally dependent.
For the \(\operatorname{Beta}(2,2)\) target,
\[ f(y)=6y(1-y),\qquad 0<y<1. \]
With a uniform proposal, \(g(y)=1\). The derivative is \(6-12y\), which is zero at \(y=1/2\); the second derivative is \(-12\) and the endpoint densities are zero. Thus
\[ M_{\min}=\max_{0<y<1}6y(1-y)=\frac{3}{2}. \]
The overall acceptance rate is \(1/M=2/3\). For a particular proposal \(y\), the decision uses \(f(y)/(Mg(y))=4y(1-y)\):
| y | u | 4y(1-y) | Decision |
|---|---|---|---|
| 0.2 | 0.5 | 0.64 | Accept; save 0.2 |
| 0.5 | 0.9 | 1.00 | Accept; save 0.5 |
| 0.9 | 0.5 | 0.36 | Reject; draw a new pair |
set.seed(8670)
y <- runif(1) # propose the value to keep
u <- runif(1) # independent decision draw
accept_prob <- 4 * y * (1 - y) # f(y) / (M * g(y))
u <= accept_prob # if TRUE, save y, not u
#> [1] TRUERepeat until \(n\) values have been saved. Choosing \(M=6\) is also valid, but the overall acceptance rate falls to \(1/6\).
set.seed(8670)
proposal <- runif(1500)
height <- 1.5 * runif(1500)
accepted <- height <= dbeta(proposal, 2, 2)
proposal_data <- data.frame(proposal, height, accepted)
ggplot(proposal_data, aes(proposal, height, color = accepted)) +
geom_point(alpha = 0.6, size = 1.2) +
stat_function(fun = function(x) dbeta(x, 2, 2),
color = "black", linewidth = 1) +
geom_hline(yintercept = 1.5, linetype = "dashed") +
scale_color_manual(values = c("grey70", "steelblue"),
labels = c("Rejected", "Accepted")) +
labs(x = "Proposal y", y = "Height", color = NULL)
The full R function in the online notes stops at exactly \(n\) accepted observations. The table compares its proposal counts with the expected cost \(nM\).
rbeta22_rejection <- function(n, M = 1.5) {
stopifnot(
length(n) == 1L, is.finite(n), n >= 0, n == floor(n),
length(M) == 1L, is.finite(M), M >= 1.5
)
draws <- numeric(n)
n_accepted <- 0L
proposals <- 0L
while (n_accepted < n) {
y <- runif(1)
u <- runif(1)
proposals <- proposals + 1L
if (u <= dbeta(y, 2, 2) / M) {
n_accepted <- n_accepted + 1L
draws[n_accepted] <- y
}
}
list(draws = draws, proposals = proposals)
}
set.seed(8670)
beta_tight <- rbeta22_rejection(2000, M = 1.5)
beta_loose <- rbeta22_rejection(2000, M = 6)
knitr::kable(booktabs = TRUE, data.frame(
M = c(1.5, 6),
proposals = c(beta_tight$proposals, beta_loose$proposals),
expected_proposals = 2000 * c(1.5, 6),
observed_acceptance = 2000 / c(beta_tight$proposals, beta_loose$proposals),
theoretical_acceptance = 1 / c(1.5, 6)
), digits = 3, col.names = c(
"M", "Proposals", "Expected proposals", "Observed rate", "Target rate"
))| M | Proposals | Expected proposals | Observed rate | Target rate |
|---|---|---|---|---|
| 1.5 | 3028 | 3000 | 0.661 | 0.667 |
| 6.0 | 11836 | 12000 | 0.169 | 0.167 |
Each proposal here uses two uniforms, one for \(Y\) and one for the decision. Thus \(n=2000\) with \(M=1.5\) uses \(3000\) proposals and \(6000\) uniform draws in expectation. For another proposal generator, its internal random-number cost could differ.
Exercise 4.2 (In Class: Check One Proposal (3 minutes)) For the same target and uniform proposal, take \(y=0.8\) and \(u=0.7\).
Solution 4.2.
The envelope must satisfy \(f(y)\leq Mg(y)\) throughout the support. A plausible histogram cannot establish that condition.
Sometimes we know only \(h(x)\) where
\[ f(x)=\frac{h(x)}{Z}, \qquad Z=\int h(x)\,dx. \]
If we can establish \(h(x)\leq Cg(x)\) everywhere, we may accept with probability \(h(Y)/\{Cg(Y)\}\) without knowing \(Z\). The acceptance rate is then \(Z/C\), not \(1/C\).
This is useful for some posterior distributions, where the normalizing constant is difficult to compute. The challenge remains finding an easy proposal and a valid global bound.
Connection to optimization. Finding the smallest bound involves maximizing \(f/g\), or its logarithm. A numerical optimizer or a finite grid may help explore this ratio, but does not by itself prove a global bound, especially on unbounded support.
Connection to importance sampling. Both methods use a proposal \(g\) and the ratio \(f/g\). Rejection sampling discards proposals to obtain draws from the target; importance sampling keeps proposals and weights their contributions when estimating an expectation.
Computational limitation. In high dimensions, a proposal that misses the target’s concentration or dependence can yield a very small acceptance rate. This helps motivate more specialized methods, including the Markov chain Monte Carlo methods encountered in other statistics courses.
Known distributional relationships are also simulation algorithms. They often reveal why a distribution appears in statistical inference.
If \(U_1,U_2\) are independent uniforms, define
\[ R=\sqrt{-2\log U_1},\qquad \Theta=2\pi U_2, \]
and
\[ Z_1=R\cos\Theta,\qquad Z_2=R\sin\Theta. \]
Then \(Z_1\) and \(Z_2\) are independent standard normal variables.
A bivariate standard normal has circularly symmetric density. Its direction is uniform, and its squared radius has a \(\chi^2_2\) distribution, equivalently an exponential distribution with rate \(1/2\). Box-Muller generates that radius and direction, then converts polar coordinates to Cartesian coordinates. The polar Jacobian contributes a factor \(r\), giving the radial density \(r e^{-r^2/2}\).
set.seed(8670)
n <- 4000
radius <- sqrt(-2 * log(runif(n)))
angle <- 2 * pi * runif(n)
z1 <- radius * cos(angle)
z2 <- radius * sin(angle)
ggplot(data.frame(z1, z2), aes(z1, z2)) +
geom_point(alpha = 0.2, color = "steelblue", size = 1) +
coord_equal() +
labs(x = expression(Z[1]), y = expression(Z[2]))
Shared ingredients do not automatically imply dependence, and an approximately zero sample correlation does not prove independence. Here independence follows from the joint density. For routine work, use rnorm(); this construction explains one way a normal generator can work.
If \(Z_1,\ldots,Z_k\) are independent \(N(0,1)\) variables,
\[ V=\sum_{j=1}^k Z_j^2\sim\chi^2_k. \]
If \(Z\sim N(0,1)\) is independent of \(V\sim\chi^2_\nu\), then
\[ T=\frac{Z}{\sqrt{V/\nu}}\sim t_\nu. \]
If \(V_1\sim\chi^2_{d_1}\) and \(V_2\sim\chi^2_{d_2}\) are independent, then
\[ \frac{V_1/d_1}{V_2/d_2}\sim F_{d_1,d_2}. \]
A \(t\) statistic divides a standardized normal estimation error by an estimated scale. When the scale estimate is unusually small, the ratio can be large. This explains the heavier tails of the \(t\) distribution and the larger critical values used when estimating an unknown normal variance.
set.seed(8670)
n <- 10000
nu <- 5
normal_numerator <- rnorm(n)
chi_square <- rowSums(matrix(rnorm(n * nu), nrow = n)^2)
t_draws <- normal_numerator / sqrt(chi_square / nu)
ggplot(data.frame(value = t_draws), aes(value)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 0.2,
fill = "lightblue", color = "white") +
stat_function(fun = function(x) dt(x, df = nu),
color = "firebrick", linewidth = 1) +
stat_function(fun = dnorm, linetype = "dashed", linewidth = 0.8) +
coord_cartesian(xlim = c(-6, 6)) +
labs(x = "Ratio", y = "Density")
The numerator and denominator must be generated independently. Reusing the same normals in both changes the distribution. In this example, \(\operatorname{Var}(T)=\nu/(\nu-2)=5/3\), so the raw \(t_5\) distribution also has greater variance than \(N(0,1)\).
If \(G_1\sim\operatorname{Gamma}(a,\lambda)\) and \(G_2\sim\operatorname{Gamma}(b,\lambda)\) are independent with the same rate, then
\[ X=\frac{G_1}{G_1+G_2}\sim\operatorname{Beta}(a,b). \]
The ratio is a random share of a positive total. This gives intuition for beta models of probabilities and proportions. The common rate cancels from the ratio.
set.seed(8670)
n <- 5000
g1 <- rgamma(n, shape = 3, rate = 1)
g2 <- rgamma(n, shape = 2, rate = 1)
beta_ratio <- g1 / (g1 + g2)
qq_data <- data.frame(
theoretical = qbeta(ppoints(n), shape1 = 3, shape2 = 2),
observed = sort(beta_ratio)
)
ggplot(qq_data, aes(theoretical, observed)) +
geom_point(alpha = 0.3, size = 0.8, color = "steelblue") +
geom_abline(intercept = 0, slope = 1, color = "firebrick") +
coord_equal() +
labs(x = "Theoretical Beta(3,2) quantile", y = "Sample quantile")
Normalizing several independent gammas with a common rate similarly produces a Dirichlet random vector. Its components sum to one, making it useful for uncertain category probabilities and mixture weights.
These identities construct new variables from familiar ones. We next distinguish two ways of combining distributions: adding contributions and selecting a population.
A sum combines contributions from all its terms:
\[ S=X_1+\cdots+X_k. \]
For independent continuous variables, the density of \(X_1+X_2\) is
\[ f_S(s)=\int f_{X_1}(x)f_{X_2}(s-x)\,dx. \]
The integral is a convolution. Simulation avoids computing this integral: generate the terms jointly according to the model and add them.
Under the independence conditions below:
size = k and prob = p.A total processing time adds durations from all stages. A count of failures before the \(k\)th success adds the failures before each successive success. The sum corresponds to an actual mechanism, not just an algebraic identity.
A finite mixture has CDF and density
\[ F_X(x)=\sum_{j=1}^K w_jF_j(x), \qquad f_X(x)=\sum_{j=1}^K w_jf_j(x), \]
where \(w_j\geq0\) and \(\sum_jw_j=1\).
Introduce an unobserved label \(C\):
\[ P(C=j)=w_j,\qquad X\mid C=j\sim F_j. \]
The law of total probability gives the mixture distribution after we average over the unknown label.
Input: component distributions \(F_1,\ldots,F_K\), weights \(w_k\geq0\) with \(\sum_k w_k=1\), and sample size \(n\).
The weights average distributions. They do not instruct us to return the numerical weighted average \(w_1X_1+\cdots+w_KX_K\).
Let \(X_1\sim N(0,1)\) and \(X_2\sim N(4,1)\) be independent. Compare:
\[ \begin{array}{ll} \text{Sum:} & S=X_1+X_2\sim N(4,2),\\ \text{Weighted average:} & A=(X_1+X_2)/2\sim N(2,1/2),\\ \text{Mixture:} & X\sim 0.5N(0,1)+0.5N(4,1). \end{array} \]
The mixture has mean 2 and variance 5. It describes observations from two subpopulations; each observation belongs to one subpopulation.
set.seed(8670)
n <- 10000
x1 <- rnorm(n, mean = 0, sd = 1)
x2 <- rnorm(n, mean = 4, sd = 1)
component <- sample.int(2, n, replace = TRUE, prob = c(0.5, 0.5))
x_mixture <- rnorm(n, mean = c(0, 4)[component], sd = 1)
mechanisms <- c("1. Sum", "2. Weighted average", "3. Mixture")
mixture_plot <- data.frame(
value = c(x1 + x2, (x1 + x2) / 2, x_mixture),
mechanism = rep(mechanisms, each = n)
)
grid <- seq(-4, 10, length.out = 500)
mixture_curves <- data.frame(
value = rep(grid, 3),
density = c(
dnorm(grid, 4, sqrt(2)),
dnorm(grid, 2, sqrt(0.5)),
0.5 * dnorm(grid, 0, 1) + 0.5 * dnorm(grid, 4, 1)
),
mechanism = rep(mechanisms, each = length(grid))
)
ggplot(mixture_plot, aes(value)) +
geom_histogram(aes(y = after_stat(density)), bins = 45,
fill = "lightblue", color = "white") +
geom_line(data = mixture_curves, aes(y = density),
color = "firebrick", linewidth = 0.9) +
facet_wrap(~ mechanism, nrow = 1) +
labs(x = "Value", y = "Density")
If component \(j\) has mean \(\mu_j\) and variance \(\sigma_j^2\), then
\[ \mathbb{E}(X)=\sum_jw_j\mu_j=\mu, \]
and the law of total variance gives
\[ \begin{aligned} \operatorname{Var}(X) &=\mathbb{E}\{\operatorname{Var}(X\mid C)\}+\operatorname{Var}\{\mathbb{E}(X\mid C)\}\\ &=\underbrace{\sum_j w_j\sigma_j^2}_{\text{within-component variation}} +\underbrace{\sum_j w_j(\mu_j-\mu)^2}_{\text{between-component variation}}. \end{aligned} \]
This decomposition explains why unobserved heterogeneity can inflate variability. A mixture of distinct normal components is generally nonnormal, but it need not be bimodal. Identical components yield the original normal distribution.
Connection to clustering and EM. Simulation starts with the label and then generates the observation. In mixture-model estimation, we see the observation and infer uncertainty about its label. This reversal is the basis for the responsibility calculations in the EM algorithm.
Exercise 4.3 (In Class: Mixture or Average? (2 minutes)) Let \(X_1\sim N(0,1)\) and \(X_2\sim N(4,1)\) be independent. T/F. Returning \((X_1+X_2)/2\) generates the same distribution as choosing one of the two components with equal probabilities. Both have mean 2; what else should you check?
Solution 4.3.
A mixture label need not take only finitely many values. Consider
\[ \Lambda\sim\operatorname{Gamma}(a,b), \qquad Y\mid\Lambda\sim\operatorname{Poisson}(\Lambda), \]
where \(b\) is a rate. Then
\[ \mathbb{E}(Y)=\frac{a}{b}, \qquad \operatorname{Var}(Y)=\frac{a}{b}+\frac{a}{b^2}. \]
The second term is variability in the underlying rate. Marginally, \(Y\) is negative binomial with size = a and prob = b / (b + 1).
set.seed(8670)
n <- 20000
a <- 2
b <- 0.5
latent_rate <- rgamma(n, shape = a, rate = b)
counts <- rpois(n, lambda = latent_rate)
knitr::kable(booktabs = TRUE, data.frame(
quantity = c("Mean", "Variance"),
simulated = c(mean(counts), var(counts)),
theoretical = c(a / b, a / b + a / b^2)
), digits = 3)| quantity | simulated | theoretical |
|---|---|---|
| Mean | 3.979 | 4 |
| Variance | 11.949 | 12 |
Customers, neighborhoods, or machines can have different event rates. Even when each unit follows a Poisson model conditional on its rate, the pooled counts may have variance greater than their mean. This is one reason negative binomial regression is useful.
We drew one independent rate for each count. If several observations share a single group-specific rate, they become dependent after that rate is integrated out. Where we place a draw in the simulation loop determines which observations share uncertainty.
The same generative logic extends to random vectors. The additional task is to reproduce how their coordinates vary together.
Suppose two measurements each have a standard normal marginal distribution. We still need to specify how they move together. Independent rnorm() calls impose independence; they cannot represent a positive correlation merely because the marginal histograms look correct.
For a normal vector, the covariance matrix supplies this information:
\[ \boldsymbol X\sim N_d(\boldsymbol\mu,\boldsymbol\Sigma). \]
If \(\boldsymbol Z\sim N_d(\boldsymbol0,I_d)\) and \(AA^\top=\boldsymbol\Sigma\), then
\[ \boldsymbol X=\boldsymbol\mu+A\boldsymbol Z. \]
Proof.
\[ \mathbb{E}(\boldsymbol X)=\boldsymbol\mu,\qquad \operatorname{Cov}(\boldsymbol X)=AI_dA^\top=\boldsymbol\Sigma. \]
The normality follows because every linear combination of \(\boldsymbol X\) is a linear combination of independent normals.
For \(n\) rows, each with target mean \(\boldsymbol\mu\) and symmetric positive-definite covariance \(\boldsymbol\Sigma\):
R <- chol(Sigma). R returns an upper triangular factor with \(R^\top R=\boldsymbol\Sigma\).\[ X=ZR+\boldsymbol1_n\boldsymbol\mu^\top. \]
For a column vector, the corresponding formula is \(X=\boldsymbol\mu+R^\top Z\). Multiplying rows by t(R) instead gives covariance \(RR^\top\), generally the wrong matrix. See the R documentation for chol().
Work one two-dimensional case before using the matrix formula. To obtain standard normal marginals with correlation \(0.6\), use
\[ \Sigma=\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}, \qquad R=\begin{pmatrix}1&0.6\\0&0.8\end{pmatrix}, \qquad R^\top R=\Sigma. \]
For independent standard normals \(Z_1,Z_2\), row multiplication gives
\[ X_1=Z_1,\qquad X_2=0.6Z_1+0.8Z_2. \]
The second variance is \(0.6^2+0.8^2=1\), and the covariance is \(0.6\operatorname{Var}(Z_1)=0.6\). If \((z_1,z_2)=(1,-1)\), the output is \((1,-0.2)\).
Sigma_small <- matrix(c(1, 0.6, 0.6, 1), nrow = 2)
R_small <- chol(Sigma_small) # upper triangular factor
z_row <- matrix(c(1, -1), nrow = 1) # one observation, two coordinates
z_row %*% R_small # one correlated observation
#> [,1] [,2]
#> [1,] 1 -0.2Exercise 4.4 (In Class: Match the Matrix to the Code (3 minutes)) Using the same \(R\), transform the row \((z_1,z_2)=(2,0)\). Multiple choice. For independent standard normal rows in Z, which gives covariance \(\Sigma\): (A) Z %*% R or (B) Z %*% t(R)?
Solution 4.4.
t(R) gives \(RR^\top=\begin{pmatrix}1.36&0.48\\0.48&0.64\end{pmatrix}\), which does not match the target.The construction has a geometric interpretation. In Figure 4.13, the same points start in a circular cloud, become an ellipse after multiplication by the matrix factor, and then move to the desired mean. Color tracks each point’s original first coordinate.
set.seed(8670)
geometry_z <- matrix(rnorm(1600), ncol = 2)
geometry_sigma <- matrix(c(1, 0.8, 0.8, 2), 2)
geometry_R <- chol(geometry_sigma)
geometry_correlated <- geometry_z %*% geometry_R
geometry_shifted <- sweep(geometry_correlated, 2, c(2, 1), FUN = "+")
stages <- c("1. Independent normals", "2. Apply the factor", "3. Add the mean")
geometry_data <- do.call(rbind, lapply(seq_along(stages), function(j) {
values <- list(geometry_z, geometry_correlated, geometry_shifted)[[j]]
data.frame(x = values[, 1], y = values[, 2],
original_z1 = geometry_z[, 1], stage = stages[j])
}))
theta <- seq(0, 2 * pi, length.out = 301)
reference_circle <- sqrt(qchisq(0.95, df = 2)) * cbind(cos(theta), sin(theta))
reference_ellipse <- reference_circle %*% geometry_R
geometry_contours <- do.call(rbind, lapply(seq_along(stages), function(j) {
values <- list(reference_circle, reference_ellipse,
sweep(reference_ellipse, 2, c(2, 1), FUN = "+"))[[j]]
data.frame(x = values[, 1], y = values[, 2], stage = stages[j])
}))
geometry_centers <- data.frame(x = c(0, 0, 2), y = c(0, 0, 1), stage = stages)
ggplot(geometry_data, aes(x, y)) +
geom_hline(yintercept = 0, color = "grey80", linewidth = 0.4) +
geom_vline(xintercept = 0, color = "grey80", linewidth = 0.4) +
geom_point(aes(color = original_z1), alpha = 0.35, size = 0.75) +
geom_path(data = geometry_contours, linewidth = 0.7) +
geom_point(data = geometry_centers, shape = 4, size = 3, stroke = 1.1) +
scale_color_gradient2(low = "#0072B2", mid = "grey80", high = "#D55E00",
midpoint = 0, guide = "none") +
facet_wrap(~ stage, nrow = 1) +
coord_equal() +
labs(x = "Coordinate 1", y = "Coordinate 2")
The reusable rmvn_chol() function in the online notes generates rows of independent normals, multiplies by chol(Sigma), and adds mu to each row. Here \(n=10000\), \(\mu=(1,2)^\top\), and \(\Sigma=\begin{pmatrix}1&0.8\\0.8&2\end{pmatrix}\). Compare the sample mean and covariance below with these targets.
rmvn_chol <- function(n, mu, Sigma) {
d <- length(mu)
stopifnot(
length(n) == 1L, is.finite(n), n >= 0, n == floor(n),
d > 0L, all(is.finite(mu)),
is.matrix(Sigma), all(dim(Sigma) == c(d, d)),
all(is.finite(Sigma)), isSymmetric(Sigma)
)
R <- chol(Sigma) # requires positive definiteness here
Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
sweep(Z %*% R, MARGIN = 2, STATS = mu, FUN = "+")
}
set.seed(8670)
mu <- c(1, 2)
Sigma <- matrix(c(1, 0.8, 0.8, 2), nrow = 2)
mvn_draws <- rmvn_chol(10000, mu, Sigma)
colMeans(mvn_draws)
#> [1] 1.017927 2.037802
cov(mvn_draws)
#> [,1] [,2]
#> [1,] 1.0040551 0.8101849
#> [2,] 0.8101849 1.9951480
Sigma
#> [,1] [,2]
#> [1,] 1.0 0.8
#> [2,] 0.8 2.0The function also returns an \(n\times d\) matrix for \(n=0\) or \(n=1\). This matters when a component in a simulated mixture receives few or no observations.
The spectral decomposition is
\[ \boldsymbol\Sigma=P\Lambda P^\top, \qquad \boldsymbol\Sigma^{1/2}=P\Lambda^{1/2}P^\top. \]
The eigenvectors specify directions, and the square roots of the eigenvalues specify how much to stretch a standard normal cloud in those directions.
This connects generation to principal component analysis. PCA identifies directions of variation in observed data; simulation constructs that variation from independent coordinates.
A covariance matrix can be positive semidefinite with some eigenvalues equal to zero. The normal vector then lies in a lower-dimensional subspace. Ordinary unpivoted Cholesky may fail, whereas an eigenvalue construction can handle that case. Tiny negative eigenvalues caused by roundoff may be set to zero after a tolerance check; substantial negative eigenvalues indicate an invalid covariance matrix and should not be silently removed.
If \(\boldsymbol\Sigma\) stays fixed across many simulation batches, compute its factor once and reuse it. Factorization costs roughly order \(d^3\), while multiplying \(n\) observations by a dense factor costs order \(nd^2\). Avoid repeating work that does not change the model.
Consider two standardized laboratory measurements with unit marginal variances. We compare independence with correlation \(\rho=0.8\) while keeping the marginals the same.
set.seed(8670)
n <- 30000
independent_pair <- rmvn_chol(n, c(0, 0), diag(2))
correlated_pair <- rmvn_chol(
n, c(0, 0), matrix(c(1, 0.8, 0.8, 1), 2)
)
joint_probability <- c(
mean(independent_pair[, 1] > 1 & independent_pair[, 2] > 1),
mean(correlated_pair[, 1] > 1 & correlated_pair[, 2] > 1)
)
joint_results <- data.frame(
model = c("Independent", "Correlation 0.8"),
probability = joint_probability,
MCSE = sqrt(joint_probability * (1 - joint_probability) / n)
)
# Plot a subset for legibility; estimates use all n observations.
dependence_plot <- rbind(
data.frame(x = independent_pair[1:2000, 1],
y = independent_pair[1:2000, 2], model = "Independent"),
data.frame(x = correlated_pair[1:2000, 1],
y = correlated_pair[1:2000, 2], model = "Correlation 0.8")
)
ggplot(dependence_plot, aes(x, y)) +
annotate("rect", xmin = 1, xmax = Inf, ymin = 1, ymax = Inf,
fill = "gold", alpha = 0.15) +
geom_point(alpha = 0.2, color = "steelblue", size = 0.8) +
geom_vline(xintercept = 1, linetype = "dashed") +
geom_hline(yintercept = 1, linetype = "dashed") +
facet_wrap(~ model) +
coord_equal() +
labs(x = "Measurement 1", y = "Measurement 2")
knitr::kable(booktabs = TRUE, joint_results, digits = 4,
col.names = c("Model", "Joint probability", "MCSE"))| Model | Joint probability | MCSE |
|---|---|---|
| Independent | 0.0261 | 0.0009 |
| Correlation 0.8 | 0.0971 | 0.0017 |
(1 - pnorm(1))^2 # exact probability under independence
#> [1] 0.02517149The simulated joint probability rises from approximately 0.026 to 0.097 under these two models. The Monte Carlo standard error (MCSE) measures simulation noise; we derive it below.
The same issue arises with simultaneous equipment failures, correlated regression errors, and repeated measurements. A study of a joint event must preserve the relevant dependence.
If \(\boldsymbol X_1,\ldots,\boldsymbol X_m\) are independent \(N_d(\boldsymbol0,\boldsymbol\Sigma)\) vectors, then
\[ W=\sum_{i=1}^m\boldsymbol X_i\boldsymbol X_i^\top \sim\operatorname{Wishart}_d(m,\boldsymbol\Sigma). \]
With observations stored in rows, this is simply crossprod(X). For a normal sample of size \(n\) with an estimated mean, its usual sample covariance satisfies
\[ (n-1)S\sim\operatorname{Wishart}_d(n-1,\boldsymbol\Sigma). \]
This explains why Wishart simulation is relevant to uncertainty in covariance estimates and multivariate inference.
set.seed(8670)
Sigma_w <- matrix(c(1, 0.4, 0.4, 2), 2)
wishart_df <- 12
W <- crossprod(rmvn_chol(wishart_df, c(0, 0), Sigma_w))
W / wishart_df # expectation equals Sigma_w, but one draw is noisy
#> [,1] [,2]
#> [1,] 0.63083420 -0.04991003
#> [2,] -0.04991003 2.89511646Each observation needs one component label for the entire vector. Drawing a different label for each coordinate changes the joint model.
set.seed(8670)
n <- 3000
component <- sample.int(2, n, replace = TRUE, prob = c(0.65, 0.35))
means <- list(c(0, 0), c(3, 2))
covariances <- list(
matrix(c(1, 0.6, 0.6, 1), 2),
matrix(c(0.6, -0.3, -0.3, 1), 2)
)
mixed_vectors <- matrix(NA_real_, nrow = n, ncol = 2)
for (j in seq_len(2)) {
rows <- which(component == j)
if (length(rows) > 0) {
mixed_vectors[rows, ] <- rmvn_chol(
length(rows), means[[j]], covariances[[j]]
)
}
}
ggplot(
data.frame(x = mixed_vectors[, 1], y = mixed_vectors[, 2],
component = factor(component)),
aes(x, y, color = component)
) +
geom_point(alpha = 0.25, size = 1) +
coord_equal() +
scale_color_manual(values = c("steelblue", "firebrick")) +
labs(x = "Variable 1", y = "Variable 2", color = "Component")
For vector means \(\boldsymbol\mu_j\) and covariances \(\boldsymbol\Sigma_j\),
\[ \operatorname{Cov}(\boldsymbol X)= \sum_jw_j\boldsymbol\Sigma_j+ \sum_jw_j(\boldsymbol\mu_j-\boldsymbol\mu) (\boldsymbol\mu_j-\boldsymbol\mu)^\top. \]
Between-group differences can create marginal association even if variables are independent within every group. This is useful when designing simulation studies for clustering or confounding.
Let \(\boldsymbol Z\sim N_d(\boldsymbol0,I_d)\). Its distribution is unchanged by rotations, so it has no preferred direction. Consequently,
\[ \boldsymbol U=\frac{\boldsymbol Z}{\|\boldsymbol Z\|} \]
is uniform on the surface of the unit sphere in \(\mathbb R^d\). The zero vector has probability zero under the ideal normal model.
Uniform direction alone puts every point at distance 1. To fill the ball uniformly in volume, independently generate a radius \(R\) satisfying
\[ P(R\leq r)=r^d,\quad 0\leq r\leq1. \]
This follows because a ball of radius \(r\) has \(r^d\) times the volume of the unit ball. Therefore,
\[ R=V^{1/d},\qquad \boldsymbol X=R\boldsymbol U, \]
where \(V\) is uniform and independent of the direction.
runif_sphere <- function(n, d) {
stopifnot(
length(n) == 1L, is.finite(n), n >= 0, n == floor(n),
length(d) == 1L, is.finite(d), d >= 1, d == floor(d)
)
Z <- matrix(rnorm(n * d), nrow = n, ncol = d)
lengths <- sqrt(rowSums(Z^2))
sweep(Z, MARGIN = 1, STATS = lengths, FUN = "/")
}We divide each row by its length directly. Constructing an \(n\times n\) diagonal matrix of inverse lengths would waste memory.
set.seed(8670)
n <- 2500
directions <- runif_sphere(n, d = 2)
disk <- sweep(directions, 1, sqrt(runif(n)), FUN = "*")
sphere_plot <- rbind(
data.frame(x = directions[, 1], y = directions[, 2],
region = "Surface: unit circle"),
data.frame(x = disk[, 1], y = disk[, 2],
region = "Interior: unit disk")
)
ggplot(sphere_plot, aes(x, y)) +
geom_point(alpha = 0.3, size = 0.7, color = "steelblue") +
facet_wrap(~ region) +
coord_equal() +
labs(x = "Coordinate 1", y = "Coordinate 2")
Random directions appear in rotationally invariant perturbations, geometric Monte Carlo methods, and randomized algorithms. In high dimensions, a uniform point in a ball is typically close to its boundary: \(P(R\leq0.9)=0.9^d\). At \(d=100\), this probability is about \(0.000027\).
Exercise 4.5 (Optional Self-Check: A Uniform Disk) T/F. A uniform angle and an independent \(R\sim U(0,1)\) give a uniform point inside the unit disk. Compare the proposed \(P(R\leq1/2)\) with the fraction of disk area inside radius \(1/2\).
Solution 4.5.
A mathematical construction tells us the target law. Numerical checks then help detect programming errors and quantify how much disagreement to expect from a finite sample.
A histogram is a useful first view, but its appearance depends on bin width and sample size. Different distributions can have the same mean, variance, or a similar-looking central region.
A convincing validation combines theory with checks that target different possible mistakes:
| Check | What it can reveal |
|---|---|
| Support and constraints | Negative waiting times, noninteger counts, proportions outside \((0,1)\) |
| Mean and variance | Incorrect location, scale, or parameterization |
| Quantiles, ECDF, and tail probabilities | Distributional errors hidden by matching moments |
| Correlations and joint events | Incorrect dependence or matrix orientation |
| Repeated runs and random-number use | Reused seeds or repeated observations |
| Acceptance rate and proposal count | A costly envelope or an implementation error |
For the rejection sampler, the ECDF should track the target CDF, and the transformed values should be approximately uniform.
ggplot(data.frame(value = beta_tight$draws), aes(value)) +
stat_ecdf(geom = "step", color = "steelblue", linewidth = 0.7) +
stat_function(fun = function(x) pbeta(x, 2, 2),
color = "firebrick", linewidth = 0.9, linetype = "dashed") +
labs(x = "Accepted value", y = "Cumulative probability")
Formal goodness-of-fit tests can supplement these checks, but passing a test is not a proof of correctness. A correct generator can occasionally be rejected, and a small sample may fail to detect a faulty one. The algorithm’s mathematical justification remains essential.
Suppose we estimate \(p=P(X>c)\) using \(B\) independent simulated observations:
\[ \hat p=\frac{1}{B}\sum_{b=1}^B\mathbf1\{X_b>c\}. \]
Each indicator is Bernoulli, so
\[ \operatorname{Var}(\hat p)=\frac{p(1-p)}{B}, \qquad \widehat{\operatorname{MCSE}}(\hat p) =\sqrt{\frac{\hat p(1-\hat p)}{B}}. \]
For an estimated mean of independent \(h(X_b)\) with finite variance, the corresponding MCSE is
\[ \widehat{\operatorname{MCSE}} \left(\frac1B\sum_{b=1}^B h(X_b)\right) =\frac{s_h}{\sqrt B}. \]
set.seed(8670)
B <- 10000
event <- rexp(B, rate = 0.7) > 3
p_hat <- mean(event)
mcse <- sqrt(p_hat * (1 - p_hat) / B)
knitr::kable(booktabs = TRUE, data.frame(
estimate = p_hat,
MCSE = mcse,
exact = exp(-0.7 * 3),
lower_MC_limit = p_hat - 1.96 * mcse,
upper_MC_limit = p_hat + 1.96 * mcse
), digits = 4,
col.names = c("Estimate", "MCSE", "Exact", "Lower MC limit", "Upper MC limit"))| Estimate | MCSE | Exact | Lower MC limit | Upper MC limit |
|---|---|---|---|---|
| 0.121 | 0.0033 | 0.1225 | 0.1146 | 0.1274 |
These limits quantify numerical uncertainty from the finite simulation. They do not quantify uncertainty about whether the exponential model or its rate describes a real population.
To halve MCSE, we generally need about four times as many independent replications. For rare events, zero simulated occurrences do not imply zero probability: the plug-in MCSE also becomes zero and is misleading. A binomial interval or a method designed for rare events is then more informative.
Sampling uncertainty concerns what would happen with another real dataset. Monte Carlo uncertainty comes from using finitely many simulated replications to study that behavior. Model uncertainty concerns whether the assumed mechanism is an adequate description of the setting.
More simulation reduces Monte Carlo uncertainty. It does not automatically reduce the other two.
| Available structure | Natural approach | Main issue to check |
|---|---|---|
| A reliable built-in generator | Use the distribution’s r function |
Parameterization and joint assumptions |
| A tractable quantile | Inverse transformation | Correct generalized inverse and endpoints |
| A tractable CDF without a convenient inverse | Numerical root finding | Bracketing, monotonicity, tolerance, and cost |
| A target density and a suitable proposal | Acceptance-rejection | Global envelope and acceptance rate |
| A known distributional identity | Transform simpler independent variables | Conditions such as independence and common rates |
| A latent-class or hierarchical model | Generate latent variables, then observations | Which observations share each latent draw |
| A normal covariance structure | Matrix factorization | Positive definiteness and factor orientation |
| Rotational or volume symmetry | Geometric construction | Surface versus volume and the correct radius |
Methods are often combined. A hierarchical model may use inverse transformation for event times, normal matrix factors for random effects, and conditional Bernoulli draws for outcomes.
The goal is a generator that is statistically correct, numerically reliable, and efficient enough for the problem. A sophisticated algorithm cannot compensate for an incorrect data-generating mechanism.
The next step is to use these generators to approximate expectations, integrals, and sampling distributions. That is the central task of Monte Carlo simulation.
The core chapter is complete. These optional examples are for independent study. Choose one to connect the generators to a full analysis; their assumptions and parameters are illustrative.
Suppose an experiment compares two independent groups, each with \(n\) observations:
\[ X_i\sim N(0,\sigma^2), \qquad Y_i\sim N(\delta,\sigma^2). \]
Here \(\delta\) is a meaningful difference we want to detect. To keep the benchmark exact, assume \(\sigma=1\) is known. The two-sided level-\(0.05\) test rejects when
\[ |Z|>z_{0.975}, \qquad Z=\frac{\bar Y-\bar X}{\sigma\sqrt{2/n}}. \]
For \(\delta=0\), the rejection probability is the Type I error rate. For \(\delta\neq0\), it is power.
We generate whole datasets, apply the proposed analysis to each dataset, and record the result. The data size \(n\) determines how informative each experiment is. The number of replications \(B\) determines how accurately we estimate its rejection probability.
set.seed(8670)
B <- 4000
sigma <- 1
design <- expand.grid(n = c(10, 20, 40, 80), delta = c(0, 0.5))
critical <- qnorm(0.975)
power_results <- do.call(rbind, lapply(seq_len(nrow(design)), function(i) {
n <- design$n[i]
delta <- design$delta[i]
control <- matrix(rnorm(B * n, 0, sigma), nrow = B)
treatment <- matrix(rnorm(B * n, delta, sigma), nrow = B)
se <- sigma * sqrt(2 / n)
statistic <- (rowMeans(treatment) - rowMeans(control)) / se
estimated <- mean(abs(statistic) > critical)
shift <- delta / se
exact <- pnorm(-critical - shift) +
pnorm(critical - shift, lower.tail = FALSE)
data.frame(
n = n, delta = delta, estimated = estimated, exact = exact,
MCSE = sqrt(estimated * (1 - estimated) / B)
)
}))
power_results$scenario <- ifelse(
power_results$delta == 0, "Null: Type I error", "Effect 0.5: power"
)
knitr::kable(booktabs = TRUE, power_results[, c("n", "delta", "estimated", "exact", "MCSE")],
digits = 4)| n | delta | estimated | exact | MCSE |
|---|---|---|---|---|
| 10 | 0.0 | 0.0503 | 0.0500 | 0.0035 |
| 20 | 0.0 | 0.0505 | 0.0500 | 0.0035 |
| 40 | 0.0 | 0.0508 | 0.0500 | 0.0035 |
| 80 | 0.0 | 0.0435 | 0.0500 | 0.0032 |
| 10 | 0.5 | 0.2020 | 0.2010 | 0.0063 |
| 20 | 0.5 | 0.3470 | 0.3526 | 0.0075 |
| 40 | 0.5 | 0.6038 | 0.6088 | 0.0077 |
| 80 | 0.5 | 0.8815 | 0.8854 | 0.0051 |
ggplot(power_results, aes(n, estimated)) +
geom_line(aes(y = exact), color = "firebrick", linewidth = 0.9) +
geom_errorbar(aes(
ymin = pmax(0, estimated - 1.96 * MCSE),
ymax = pmin(1, estimated + 1.96 * MCSE)
), width = 3, color = "steelblue") +
geom_point(color = "steelblue", size = 2) +
facet_wrap(~ scenario) +
scale_x_continuous(breaks = c(10, 20, 40, 80)) +
coord_cartesian(ylim = c(0, 1)) +
labs(x = "Sample size per group", y = "Rejection probability")
The exact formula gives a benchmark for our simulation. When the model becomes more realistic, such as unequal variances, clustered observations, or missing outcomes, simulation may remain feasible even when a simple power formula does not.
If \(\sigma\) is estimated, replace the known-variance test with the actual planned procedure, such as a two-sample \(t\) test. Simulating one analysis and reporting power for another answers the wrong question.
Suppose a process produces independent binary defect indicators conditional on an unknown defect probability \(p\). Start with a \(\operatorname{Beta}(1,1)\) prior and observe 3 defects among 18 items. Beta-binomial conjugacy gives
\[ p\mid\text{data}\sim\operatorname{Beta}(4,16). \]
To predict the number of defects \(Y_{\mathrm{new}}\) in a future batch of \(m=50\) items, generate
\[ p^{(b)}\sim\operatorname{Beta}(4,16), \qquad Y_{\mathrm{new}}^{(b)}\mid p^{(b)} \sim\operatorname{Binomial}(50,p^{(b)}). \]
Each simulated future batch uses one shared draw of its possible process rate. This integrates uncertainty about \(p\) into the prediction.
The plug-in alternative uses only \(\bar p=4/20=0.2\):
\[ Y_{\mathrm{plug-in}}\sim\operatorname{Binomial}(50,0.2). \]
Both have mean 10, but the posterior predictive variance is larger:
\[ \begin{aligned} \operatorname{Var}(Y_{\mathrm{new}}\mid\text{data}) &=m\bar p(1-\bar p)\frac{a+b+m}{a+b+1}\\ &=50(0.2)(0.8)\frac{70}{21} \approx26.67, \end{aligned} \]
compared with the plug-in variance \(50(0.2)(0.8)=8\).
set.seed(8670)
B <- 20000
a_post <- 4
b_post <- 16
batch_size <- 50
p_mean <- a_post / (a_post + b_post)
p_draw <- rbeta(B, shape1 = a_post, shape2 = b_post)
predictive_defects <- rbinom(B, size = batch_size, prob = p_draw)
plugin_defects <- rbinom(B, size = batch_size, prob = p_mean)
k <- 0:batch_size
# Exact beta-binomial probabilities, computed on the log scale.
predictive_pmf <- exp(
lchoose(batch_size, k) +
lbeta(k + a_post, batch_size - k + b_post) -
lbeta(a_post, b_post)
)
predictive_comparison <- data.frame(
model = c("Posterior predictive", "Plug-in"),
simulated_variance = c(var(predictive_defects), var(plugin_defects)),
theoretical_variance = c(
batch_size * p_mean * (1 - p_mean) *
(a_post + b_post + batch_size) / (a_post + b_post + 1),
batch_size * p_mean * (1 - p_mean)
),
simulated_prob_at_least_15 = c(
mean(predictive_defects >= 15), mean(plugin_defects >= 15)
),
exact_prob_at_least_15 = c(
sum(predictive_pmf[k >= 15]),
pbinom(14, size = batch_size, prob = p_mean, lower.tail = FALSE)
)
)
knitr::kable(booktabs = TRUE, predictive_comparison, digits = 3, col.names = c(
"Model", "Sim. Var", "Exact Var", "Sim. tail", "Exact tail"
))| Model | Sim. Var | Exact Var | Sim. tail | Exact tail |
|---|---|---|---|---|
| Posterior predictive | 27.065 | 26.667 | 0.194 | 0.187 |
| Plug-in | 8.004 | 8.000 | 0.061 | 0.061 |
The tail columns report \(P(Y\geq15)\).
predictive_plot <- rbind(
data.frame(
defects = k,
proportion = tabulate(predictive_defects + 1, nbins = batch_size + 1) / B,
exact = predictive_pmf,
model = "Posterior predictive"
),
data.frame(
defects = k,
proportion = tabulate(plugin_defects + 1, nbins = batch_size + 1) / B,
exact = dbinom(k, size = batch_size, prob = p_mean),
model = "Plug-in"
)
)
ggplot(predictive_plot, aes(defects, proportion)) +
geom_col(fill = "lightblue", width = 0.85) +
geom_line(aes(y = exact), color = "firebrick", linewidth = 0.8) +
geom_vline(xintercept = 14.5, linetype = "dashed") +
facet_wrap(~ model) +
coord_cartesian(xlim = c(0, 35)) +
labs(x = "Defects in the next 50 items", y = "Probability")
A decision triggered by a large defect count depends on the predictive tail, not just the expected count. Hierarchical simulation carries parameter uncertainty through to that tail.
If we drew a fresh independent \(p\) for each individual item instead of one shared \(p\) for the batch, the marginal total would be binomial with probability \(\mathbb{E}(p)\). That would represent a different mechanism and lose the shared uncertainty in this posterior prediction.
Consider a single server, first-come-first-served, initially idle. Let interarrival times be independent exponentials with rate \(\lambda\), and service times be independent \(\operatorname{Gamma}(2,1)\) variables, with mean 2 minutes. Assume the arrival and service processes are independent.
If \(A_i\) is the time between customers \(i-1\) and \(i\), \(S_i\) is customer \(i\)’s service time, and \(W_i\) is the wait before service, then
\[ W_1=0,\qquad W_i=\max\{0,W_{i-1}+S_{i-1}-A_i\}. \]
This recursion turns independent primitive random inputs into dependent waiting times. A long service can delay several later customers.
We compare arrival rates \(0.35\) and \(0.45\) per minute. Their utilization values are \(\lambda\mathbb{E}(S)=0.7\) and \(0.9\). Both are below one, but the system closer to capacity has less room to absorb bursts.
simulate_queue <- function(arrival_rate, n_customers = 200) {
interarrival <- rexp(n_customers, rate = arrival_rate)
service <- rgamma(n_customers, shape = 2, rate = 1)
wait <- numeric(n_customers)
if (n_customers > 1) {
for (i in 2:n_customers) {
wait[i] <- max(0, wait[i - 1] + service[i - 1] - interarrival[i])
}
}
wait
}
set.seed(8670)
B_days <- 800
rates <- c(0.35, 0.45)
queue_results <- do.call(rbind, lapply(rates, function(rate) {
# Each replication is an independent run of 200 customers, starting empty.
fraction_delayed <- replicate(
B_days, mean(simulate_queue(rate) > 5)
)
data.frame(
arrival_rate = rate,
utilization = 2 * rate,
fraction_waiting_over_5 = mean(fraction_delayed),
MCSE = sd(fraction_delayed) / sqrt(B_days)
)
}))
knitr::kable(booktabs = TRUE, queue_results, digits = 4, col.names = c(
"Arrival rate", "Utilization", "Fraction waiting > 5", "MCSE"
))| Arrival rate | Utilization | Fraction waiting > 5 | MCSE |
|---|---|---|---|
| 0.35 | 0.7 | 0.2385 | 0.0037 |
| 0.45 | 0.9 | 0.5642 | 0.0066 |
set.seed(8670)
queue_paths <- do.call(rbind, lapply(rates, function(rate) {
data.frame(
customer = 1:200, wait = simulate_queue(rate),
scenario = paste("Arrival rate =", rate)
)
}))
ggplot(queue_paths, aes(customer, wait)) +
geom_line(color = "steelblue", linewidth = 0.6) +
geom_hline(yintercept = 5, linetype = "dashed", color = "firebrick") +
facet_wrap(~ scenario) +
labs(x = "Customer index", y = "Wait before service (minutes)")
The estimand is the expected fraction of the first 200 customers who wait more than five minutes in a system starting empty. It is not automatically the long-run stationary probability. Studying steady-state performance requires addressing the initial transient and dependence.
We compute MCSE from independent run-level fractions. Pooling all customers and treating their wait indicators as independent Bernoulli trials would ignore the dependence created by the queue.
This example combines distribution choice, recursive computation, and statistical reasoning about the unit of replication.
Try the concept checks first (about 5 minutes), then the short calculations (about 10 minutes). Explain your choice before opening the solution.
Exercise 4.6 (Concept Check: Choose or Mark True/False)
Solution 4.6.
| Item | Answer | Reason |
|---|---|---|
| 1 | B | Apply \(F^{-1}\) to convert the probability level into an observation. |
| 2 | F | The rate is \(1/M\) for a normalized target and proposal. |
| 3 | F | \(Z\) is unnecessary for the decision; the overall rate is \(Z/C\). |
| 4 | A | Divide by the square root of the independent chi-square variable per degree of freedom. |
| 5 | T | Both coordinates belong to the chosen component; separate labels change the joint model. |
| 6 | F | Marginals do not determine dependence; also check covariance and joint events. |
Exercise 4.7 (Practice: Construct, Calculate, Check)
Z contain independent standard normals. The target has mean mu and positive definite covariance Sigma. Correct the second line and explain which covariance check detects the original error.R <- chol(Sigma)
X <- sweep(Z %*% t(R), 2, mu, FUN = "+")Solution 4.7.
Inverse transformation (follow Figure 4.1).
Rejection sampling (follow Figure 4.5).
Repair the R code (follow Figure 4.12).
chol(Sigma) gives \(R^\top R=\Sigma\).X <- sweep(Z %*% R, 2, mu, FUN = "+"). This gives covariance \(R^\top R\); the original gives \(RR^\top\).cov(X) with Sigma, allowing for Monte Carlo variation. Adding mu changes the mean, not covariance.This chapter adapts and expands the 2025 STAT 8670 chapter on generating random variables and its Quarto source in the master branch.
For software details, consult the official R documentation linked at the relevant examples. Useful local help pages include ?RNGkind, ?Distributions, ?qgeom, ?chol, and ?rWishart. The course’s original reading is Rizzo, M. L. (2007), Statistical Computing with R, Chapman & Hall/CRC.