In Chapter 3: Optimization, 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?) 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.
A normalized density has total area 1. An unnormalized target has the same shape but need not have area 1. Either form works for acceptance-rejection sampling with a valid envelope. Further details are in Supplementary (not on the exam).
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)
u <- runif(1500)
v <- 1.5 * u
accepted <- v <= dbeta(proposal, 2, 2)
proposal_data <- data.frame(proposal, v, accepted)
accepted_data <- data.frame(value = proposal[accepted])
target_data <- data.frame(value = seq(0, 1, length.out = 501))
target_data$density <- dbeta(target_data$value, 2, 2)
left_plot <- ggplot(proposal_data, aes(proposal, v)) +
geom_point(aes(color = ifelse(accepted, "Accepted", "Rejected")),
alpha = 0.6, size = 1.2) +
geom_hline(yintercept = 1.5, linetype = "dashed", color = "grey40") +
geom_line(data = target_data, aes(value, density),
color = "black", linewidth = 1, linetype = "longdash") +
scale_color_manual(values = c("Rejected" = "grey70", "Accepted" = "steelblue"),
breaks = c("Rejected", "Accepted")) +
scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
scale_y_continuous(limits = c(0, 1.75), breaks = seq(0, 1.6, by = 0.4),
expand = expansion(mult = c(0, 0))) +
labs(title = "Acceptance-rejection", x = "Proposal y",
y = "Auxiliary ordinate v\n(density scale)", color = NULL) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
theme(legend.position = "bottom", legend.text = element_text(size = 11))
right_plot <- ggplot() +
geom_histogram(data = accepted_data, aes(value, after_stat(density)),
binwidth = 0.05, boundary = 0,
fill = "steelblue", color = "white", alpha = 0.25) +
geom_density(data = accepted_data,
aes(value, color = "Accepted-sample density"),
bounds = c(0, 1), linewidth = 1, key_glyph = "path") +
geom_line(data = target_data, aes(value, density, color = "Target Beta(2,2)"),
linewidth = 1, linetype = "longdash") +
scale_color_manual(values = c("Accepted-sample density" = "steelblue",
"Target Beta(2,2)" = "black"),
breaks = c("Accepted-sample density", "Target Beta(2,2)")) +
scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25)) +
scale_y_continuous(limits = c(0, 1.75), breaks = seq(0, 1.6, by = 0.4),
expand = expansion(mult = c(0, 0))) +
labs(title = "Accepted sample vs. target", x = "Accepted y", y = "Density", color = NULL) +
guides(color = guide_legend(nrow = 2, byrow = TRUE)) +
theme(legend.position = "bottom", legend.text = element_text(size = 11))
# Use separate plots so the two vertical axes describe different quantities.
grid::grid.newpage()
grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2)))
print(left_plot, newpage = FALSE,
vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1))
print(right_plot, newpage = FALSE,
vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2))
grid::popViewport()
Save the horizontal coordinates \(Y\) of the accepted blue points. The vertical heights are used only for the accept/reject decision.
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) 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.
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),
aes(color = "t distribution (df = 5)",
linetype = "t distribution (df = 5)"), linewidth = 1) +
stat_function(fun = dnorm,
aes(color = "Standard normal N(0, 1)",
linetype = "Standard normal N(0, 1)"), linewidth = 0.8) +
scale_color_manual(name = "Density curve",
breaks = c("t distribution (df = 5)", "Standard normal N(0, 1)"),
values = c("t distribution (df = 5)" = "firebrick",
"Standard normal N(0, 1)" = "black")) +
scale_linetype_manual(name = "Density curve",
breaks = c("t distribution (df = 5)", "Standard normal N(0, 1)"),
values = c("t distribution (df = 5)" = "solid",
"Standard normal N(0, 1)" = "dashed")) +
coord_cartesian(xlim = c(-6, 6)) +
labs(x = "Ratio", y = "Density") +
theme(legend.position = "bottom", legend.text = element_text(size = 11))
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")
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.
The first two panels of Figure 4.11 show the component densities separately. In the third panel, multiply each density by its selection probability and add the two curves. The mixture retains observations near both component centers; it does not move each observation to their average.
mixture_panels <- c("1. Component 1: N(0, 1)",
"2. Component 2: N(4, 1)",
"3. Mixture: 0.5 f1 + 0.5 f2")
mixture_grid <- seq(-4, 8, length.out = 601)
mixture_f1 <- dnorm(mixture_grid, 0, 1)
mixture_f2 <- dnorm(mixture_grid, 4, 1)
mixture_densities <- data.frame(
x = rep(mixture_grid, 3),
density = c(mixture_f1, mixture_f2, 0.5 * mixture_f1 + 0.5 * mixture_f2),
panel = factor(rep(mixture_panels, each = length(mixture_grid)),
levels = mixture_panels),
curve = rep(c("Component 1", "Component 2", "Mixture"),
each = length(mixture_grid))
)
mixture_contributions <- data.frame(
x = rep(mixture_grid, 2), density = c(0.5 * mixture_f1, 0.5 * mixture_f2),
panel = factor(mixture_panels[3], levels = mixture_panels),
curve = rep(c("Component 1", "Component 2"), each = length(mixture_grid))
)
mixture_colors <- c("Component 1" = "#0072B2", "Component 2" = "#D55E00",
"Mixture" = "#303030")
ggplot(mixture_densities, aes(x, density)) +
geom_area(aes(fill = curve), alpha = 0.12) +
geom_line(data = mixture_contributions, aes(color = curve),
linetype = "dashed", linewidth = 0.8) +
geom_line(aes(color = curve), linewidth = 1) +
scale_color_manual(values = mixture_colors, guide = "none") +
scale_fill_manual(values = mixture_colors, guide = "none") +
scale_x_continuous(breaks = c(-4, 0, 4, 8)) +
scale_y_continuous(limits = c(0, 0.43), breaks = c(0, 0.2, 0.4)) +
facet_wrap(~ panel, nrow = 1) +
labs(x = "Value", y = "Density") +
theme(strip.text = element_text(size = 11, face = "bold"),
panel.grid.minor = element_blank())
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?) 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.
We now extend simulation to random vectors, where we must also reproduce dependence between coordinates.
Correct marginal distributions are not enough: we must also reproduce how the coordinates vary together. To generate \(\boldsymbol X\sim N_d(\boldsymbol\mu,\boldsymbol\Sigma)\), start with independent standard normals and introduce the required covariance through a matrix transformation.
If \(\boldsymbol Z\sim N_d(\boldsymbol0,I_d)\) and \(AA^\top=\boldsymbol\Sigma\), then
\[ \boldsymbol X=\boldsymbol\mu+A\boldsymbol Z \sim N_d(\boldsymbol\mu,\boldsymbol\Sigma). \]
The transformation preserves normality, with \(E(\boldsymbol X)=\boldsymbol\mu\) and \(\operatorname{Cov}(\boldsymbol X)=AA^\top=\boldsymbol\Sigma\).
For a symmetric positive-definite Sigma, R’s R <- chol(Sigma) returns an upper triangular factor satisfying \(R^\top R=\Sigma\). When observations are stored in rows:
Z of independent standard normals.Z %*% R to introduce covariance \(\Sigma\).mu to each row with sweep(..., 2, mu, "+").Using Z %*% t(R) gives \(RR^\top\) instead, generally the wrong covariance. See the R documentation for chol().
Suppose the target mean and covariance are
\[ \boldsymbol\mu=\begin{pmatrix}1\\2\end{pmatrix},\qquad \Sigma=\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix},\qquad R=\begin{pmatrix}1&0.6\\0&0.8\end{pmatrix}. \]
Here the transformation is simply
\[ X_1=1+Z_1,\qquad X_2=2+0.6Z_1+0.8Z_2. \]
Both variances are 1 because \(0.6^2+0.8^2=1\), and the covariance is \(0.6\). Thus the correlation is also \(0.6\).
set.seed(8670)
n <- 10000
mu <- c(1, 2)
Sigma <- matrix(c(1, 0.6, 0.6, 1), nrow = 2)
R <- chol(Sigma)
Z <- matrix(rnorm(n * 2), nrow = n, ncol = 2)
X <- sweep(Z %*% R, 2, mu, FUN = "+")
colMeans(X) # should be close to c(1, 2)
#> [1] 1.017927 2.026850
cov(X) # should be close to Sigma
#> [,1] [,2]
#> [1,] 1.0040551 0.6071944
#> [2,] 0.6071944 0.9984428The sample checks should be close to their targets, with some simulation error. In Figure 4.12, multiplication introduces correlation and adding the mean shifts the cloud.
geometry_z <- Z[seq_len(800), , drop = FALSE]
geometry_correlated <- geometry_z %*% R
geometry_shifted <- X[seq_len(800), , drop = FALSE]
stages <- c("1. Independent normals", "2. Introduce correlation", "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 %*% R
geometry_contours <- do.call(rbind, lapply(seq_along(stages), function(j) {
values <- list(reference_circle, reference_ellipse,
sweep(reference_ellipse, 2, mu, FUN = "+"))[[j]]
data.frame(x = values[, 1], y = values[, 2], stage = stages[j])
}))
geometry_centers <- data.frame(x = c(0, 0, mu[1]), y = c(0, 0, mu[2]), 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")
Exercise 4.4 (In Class: Match the Matrix to the Code) Using this example’s \(R\) and \(\mu\), transform the row \((z_1,z_2)=(2,0)\) and add the mean. Multiple choice. Which multiplication gives covariance \(\Sigma\): (A) Z %*% R or (B) Z %*% t(R)?
Solution 4.4. The transformed row is \((2,0)R+(1,2)=(3,3.2)\). Choose A, since the covariance of each transformed row is \(R^\top R=\Sigma\).
The same two-step rule works for vectors. For a two-component mixture, choose \(C=1\) with probability \(w\) and \(C=2\) with probability \(1-w\), then generate the whole vector from the selected multivariate normal distribution:
\[ \boldsymbol X\mid C=j\sim N_d(\boldsymbol\mu_j,\Sigma_j), \qquad f(\boldsymbol x)=w f_1(\boldsymbol x)+(1-w)f_2(\boldsymbol x). \]
Here \(f_1\) and \(f_2\) are joint densities. Use one label for all coordinates of an observation and a fresh label for the next observation. Choosing a separate label for each coordinate generally gives a different joint distribution. A mixture of multivariate normals need not itself be multivariate normal.
Consider equal weights and the two components used in last year’s notes:
\[ \boldsymbol\mu_1=\begin{pmatrix}0\\0\end{pmatrix},\quad \Sigma_1=\begin{pmatrix}1&0\\0&1\end{pmatrix},\qquad \boldsymbol\mu_2=\begin{pmatrix}3\\2\end{pmatrix},\quad \Sigma_2=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}. \]
Component 1 has independent coordinates and circular density contours. Component 2 has correlation 0.5 and tilted contours. To simulate the mixture, draw a label for each observation, then use that component’s mean and Cholesky factor to generate the whole row.
set.seed(8670)
n <- 1000
component <- sample.int(2, n, replace = TRUE, prob = c(0.5, 0.5))
means <- list(c(0, 0), c(3, 2))
covariances <- list(diag(2), matrix(c(1, 0.5, 0.5, 1), 2, 2))
X_mix <- matrix(NA_real_, nrow = n, ncol = 2)
for (j in 1:2) {
rows <- which(component == j)
Z <- matrix(rnorm(2 * length(rows)), ncol = 2)
X_mix[rows, ] <- sweep(Z %*% chol(covariances[[j]]), 2, means[[j]], "+")
}Each row of X_mix is one bivariate observation. In Figure 4.13, blue circles come from component 1 and orange triangles come from component 2. The first two panels separate the observations by label. The right panel combines those same points, keeping their coordinates and colors unchanged.
The gray contours show the component densities and the combined joint density of \((X_1,X_2)\). At each point, the mixture density is half of the first component’s density plus half of the second’s. Colors identify the generating component, even where the two clouds overlap.
joint_xlim <- range(c(-3.5, 6.5, X_mix[, 1]))
joint_ylim <- range(c(-3.5, 5.5, X_mix[, 2]))
joint_grid <- expand.grid(x1 = seq(joint_xlim[1], joint_xlim[2], length.out = 201),
x2 = seq(joint_ylim[1], joint_ylim[2], length.out = 181))
# Both component variances are 1; rho controls their correlation.
bivariate_density <- function(x1, x2, mu, rho) {
z1 <- x1 - mu[1]
z2 <- x2 - mu[2]
exp(-(z1^2 - 2 * rho * z1 * z2 + z2^2) / (2 * (1 - rho^2))) /
(2 * pi * sqrt(1 - rho^2))
}
joint_f1 <- with(joint_grid, bivariate_density(x1, x2, c(0, 0), 0))
joint_f2 <- with(joint_grid, bivariate_density(x1, x2, c(3, 2), 0.5))
joint_panels <- c("1. Component 1: f1", "2. Component 2: f2",
"3. Joint mixture: 0.5 f1 + 0.5 f2")
joint_densities <- do.call(rbind, lapply(seq_along(joint_panels), function(j) {
data.frame(joint_grid,
density = list(joint_f1, joint_f2, 0.5 * joint_f1 + 0.5 * joint_f2)[[j]],
panel = factor(joint_panels[j], levels = joint_panels))
}))
# Reuse each observation in its component panel and in the mixture panel.
joint_observations <- data.frame(
id = seq_len(n), x1 = X_mix[, 1], x2 = X_mix[, 2],
source = factor(component, levels = 1:2,
labels = c("Component 1", "Component 2"))
)
joint_scatter <- rbind(
transform(joint_observations, panel = joint_panels[component]),
transform(joint_observations, panel = joint_panels[3])
)
joint_scatter$panel <- factor(joint_scatter$panel, levels = joint_panels)
ggplot(joint_densities, aes(x1, x2)) +
geom_raster(aes(fill = density)) +
geom_contour(aes(z = density), breaks = c(0.01, 0.03, 0.06, 0.10, 0.15),
color = "#777777", linewidth = 0.4) +
geom_point(data = joint_scatter, aes(color = source, shape = source),
size = 0.95, alpha = 0.55) +
scale_fill_gradient(low = "#FFFFFF", high = "#E0E0E0",
limits = c(0, 0.19), guide = "none") +
scale_color_manual(values = c("Component 1" = "#0072B2",
"Component 2" = "#D55E00"), name = NULL) +
scale_shape_manual(values = c("Component 1" = 16, "Component 2" = 17),
name = NULL) +
guides(color = guide_legend(override.aes = list(alpha = 1, size = 2.8))) +
scale_x_continuous(breaks = c(-3, 0, 3, 6), expand = c(0, 0)) +
scale_y_continuous(breaks = c(-3, 0, 3, 5), expand = c(0, 0)) +
facet_wrap(~ panel, nrow = 1) +
coord_equal() +
labs(x = expression(X[1]), y = expression(X[2])) +
theme(strip.text = element_text(size = 10, face = "bold"),
legend.position = "bottom",
panel.grid = element_blank())
Start with a few simple checks against the intended distribution:
| Check | What to compare |
|---|---|
| Possible values | Are all draws in the correct range? Counts should be integers; exponential waiting times should be nonnegative. |
| Mean and variance | Compare the sample summaries with the theoretical mean and variance. |
| Distributional shape | Compare a histogram with the target density. |
| Multivariate dependence | Inspect a scatterplot and compare the sample covariance or correlation with the target. |
For example, an exponential distribution with rate 2 has mean \(1/2\) and variance \(1/4\):
set.seed(8670)
x <- rexp(1000, rate = 2)
range(x) # all values should be nonnegative
mean(x) # should be close to 1/2
var(x) # should be close to 1/4
hist(x, probability = TRUE)
curve(dexp(x, rate = 2), add = TRUE, col = "red")Small differences are expected in a random sample. These checks can reveal errors, but matching a few summaries does not prove that a generator is correct. If a check fails, review the formula, parameters, and code; generating more observations does not repair an incorrect algorithm.
Try the concept checks first, then the short calculations. Explain your choice before opening the solution.
Exercise 4.5 (Concept Check: Choose or Mark True/False)
Solution 4.5.
| 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 | Count accepted values only; continue sampling after a rejection. |
| 4 | A | Divide by the square root of the independent chi-square variable per degree of freedom. |
| 5 | T | Independently select a component for each observation, then draw from that component. |
| 6 | F | Marginals do not determine dependence; also inspect a scatterplot and compare covariance or correlation. |
Exercise 4.6 (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.6.
Inverse transformation (follow Figure 4.1).
Rejection sampling (follow Figure 4.5).
Repair the R code (use the Cholesky construction above).
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 answers one question: How can we turn a probability model into simulated data?
The content here won’t be on the exam
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.
The core chapter is complete. These 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.
These optional details will not be on the exam. They give the derivations and plots behind the brief comparison in the main text.
Each point in the left panel of Figure 4.6 is \((Y,V)\). Generate a proposal \(Y\sim\operatorname{Uniform}(0,1)\) and an independent decision variable \(U\sim\operatorname{Uniform}(0,1)\), then plot \(V=Mg(Y)U=1.5U\). Thus the plotted ordinate ranges from 0 to 1.5, while the original decision variable \(U\) ranges from 0 to 1. The two acceptance rules are equivalent:
\[ V\leq f(Y) \quad\Longleftrightarrow\quad U\leq\frac{f(Y)}{Mg(Y)}=4Y(1-Y). \]
Blue points satisfy this rule. Save their horizontal coordinates \(Y\); discard \(U\) and \(V\). A point’s ordinate is an auxiliary random value used for the decision, not \(f(Y)\), an estimated density, or an acceptance probability. The left axis therefore says Auxiliary ordinate v (density scale). The right panel estimates the density of the retained \(Y\) values, so its axis says Density. The horizontal envelope at 1.5 is a bound on \(f\), not a probability density with area 1.
An unnormalized target is a nonnegative function \(q(y)\) whose relative heights describe the desired distribution. Its total area is \(Z=\int q(y)\,dy\), where \(0<Z<\infty\). The normalized density is \(f(y)=q(y)/Z\), so \(\int f(y)\,dy=1\). Dividing by the area preserves the shape; it does not set the maximum to 1.
For the Beta(2,2) example, choose \(q(y)=y(1-y)\) on \([0,1]\) and zero elsewhere. Then
\[ Z=\int_0^1 y(1-y)\,dy =\left[\frac{y^2}{2}-\frac{y^3}{3}\right]_0^1 =\frac16, \qquad f(y)=\frac{q(y)}{1/6}=6y(1-y). \]
| Quantity | Unnormalized target \(q\) | Normalized density \(f=q/Z\) |
|---|---|---|
| Function on \([0,1]\) | \(y(1-y)\) | \(6y(1-y)\) |
| Total area | \(1/6\) | \(1\) |
| Maximum at \(y=1/2\) | \(1/4\) | \(3/2=1.5\) |
Why use the normalized version in Figure 4.6? The right-hand histogram uses bar height \(\text{count}/(n_{\text{accepted}}\times\text{bin width})\), so its bar areas sum to 1; the smooth density estimate also has area 1. Overlaying \(f\), which also has area 1, makes a direct comparison possible. Overlaying \(q\) would put the target six times too low even when the accepted sample is correct.
Probabilities are areas over intervals: \(P(a\leq Y\leq b)=\int_a^b f(y)\,dy\). For example, a histogram density of \(1.6\) over a bin of width \(0.05\) represents relative frequency \(1.6\times0.05=0.08\). The target is already normalized even though its peak is 1.5.
Normalization is needed for a density comparison, but the sampling algorithm need not know \(Z\). If the proposal density \(g\) is normalized and a known constant \(C\) satisfies \(q(y)\leq Cg(y)\) everywhere, accept with probability \(q(y)/(Cg(y))\). Then
\[ P(\text{accept})=\frac{Z}{C}, \qquad f_{Y\mid\text{accept}}(y) =\frac{g(y)\,q(y)/(Cg(y))}{Z/C} =\frac{q(y)}{Z}=f(y). \]
Here \(g(y)=1\), \(q(y)=y(1-y)\) and \(C=1/4\), giving acceptance probability \(4y(1-y)\) and overall rate \((1/6)/(1/4)=2/3\). These are exactly the same as using normalized \(f\) with \(M=3/2\). The corresponding unnormalized plot would use ordinate \(W=CU=U/4\) and compare it with \(q(Y)\); both the target and envelope are one-sixth of the scales used in Figure 4.6. Keep the target and envelope on the same scale. The formula \(1/M\) for the overall acceptance rate assumes a normalized target; with \(q\), the rate is \(Z/C\).
\[ \frac{f(Y)}{3/2}=\frac{q(Y)}{1/4}=4Y(1-Y). \]
Keeping \(3/2\) as the envelope for \(q\) is still valid here, but less efficient. Do not compare \(q\) directly with a density histogram: their vertical scales do not match.
Figure 4.17 uses the same proposal and decision pairs \((Y,U)\) in both panels. The right panel divides the target, envelope, and plotted ordinate by 6. Consequently, every point keeps the same accept/reject classification:
\[ U\leq\frac{f(Y)}{3/2} =\frac{q(Y)}{1/4} =4Y(1-Y). \]
set.seed(8670)
comparison_pairs <- data.frame(y = runif(160), u = runif(160))
comparison_pairs$decision <- ifelse(
comparison_pairs$u <= 4 * comparison_pairs$y * (1 - comparison_pairs$y),
"Accepted", "Rejected"
)
comparison_examples <- data.frame(
y = c(0.2, 0.8), u = c(0.4, 0.8),
decision = c("Accepted", "Rejected"),
label = c("A: accept", "B: reject")
)
comparison_grid <- data.frame(y = seq(0, 1, length.out = 501))
comparison_colors <- c("Accepted" = "steelblue", "Rejected" = "grey55")
comparison_panel <- function(envelope, target_multiplier, title,
target_label, envelope_label, axis_label,
acceptance_label, y_breaks) {
points <- transform(comparison_pairs, ordinate = envelope * u)
examples <- transform(comparison_examples, ordinate = envelope * u)
curve <- transform(comparison_grid,
target = target_multiplier * y * (1 - y))
ggplot() +
geom_ribbon(data = curve, aes(y, ymin = 0, ymax = target),
fill = "steelblue", alpha = 0.12) +
geom_point(data = points, aes(y, ordinate, color = decision),
alpha = 0.55, size = 1.5) +
geom_line(data = curve, aes(y, target), color = "black", linewidth = 1) +
geom_hline(yintercept = envelope, color = "#B45F06",
linetype = "dashed", linewidth = 0.9) +
annotate("label", x = 0.5, y = 0.28 * envelope, label = target_label,
size = 3.6, color = "black", fill = "white", linewidth = 0) +
annotate("label", x = 0.5, y = 1.13 * envelope, label = envelope_label,
size = 3.7, color = "#9A4E00", fill = "white", linewidth = 0,
lineheight = 1.1) +
geom_point(data = examples, aes(y, ordinate, fill = decision),
shape = 21, color = "black", size = 3.3, stroke = 0.8,
show.legend = FALSE) +
geom_label(data = examples,
aes(y, ordinate + 0.09 * envelope, label = label),
size = 3.4, fill = "white", linewidth = 0.2) +
scale_color_manual(values = comparison_colors, breaks = names(comparison_colors)) +
scale_fill_manual(values = comparison_colors, guide = "none") +
scale_x_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.25),
expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(limits = c(0, 1.25 * envelope), breaks = y_breaks,
expand = expansion(mult = c(0, 0))) +
labs(title = title, x = "Same proposal y", y = axis_label,
color = NULL, caption = acceptance_label) +
theme(legend.position = "bottom", legend.text = element_text(size = 11),
plot.title = element_text(size = 16),
plot.caption = element_text(size = 11, hjust = 0.5),
plot.margin = margin(8, 10, 4, 8))
}
normalized_comparison <- comparison_panel(
envelope = 1.5, target_multiplier = 6, title = "Normalized target",
target_label = "Target f(y) = 6y(1-y)\nTotal area = 1",
envelope_label = "Envelope Mg(y) = 3/2\nTotal area = 3/2",
axis_label = "Auxiliary ordinate V = 1.5U\n(normalized density scale)",
acceptance_label = "Acceptance probability = 1 / (3/2) = 2/3",
y_breaks = c(0, 0.3, 0.6, 0.9, 1.2, 1.5)
)
unnormalized_comparison <- comparison_panel(
envelope = 0.25, target_multiplier = 1, title = "Unnormalized target",
target_label = "Target q(y) = y(1-y)\nTotal area = 1/6",
envelope_label = "Envelope Cg(y) = 1/4\nTotal area = 1/4",
axis_label = "Auxiliary ordinate W = U/4\n(unnormalized target scale)",
acceptance_label = "Acceptance probability = (1/6) / (1/4) = 2/3",
y_breaks = seq(0, 0.25, by = 0.05)
)
# Both scales classify every shared (Y,U) pair identically.
stopifnot(identical(
1.5 * comparison_pairs$u <= dbeta(comparison_pairs$y, 2, 2),
0.25 * comparison_pairs$u <= comparison_pairs$y * (1 - comparison_pairs$y)
))
grid::grid.newpage()
grid::pushViewport(grid::viewport(layout = grid::grid.layout(
2, 2, heights = grid::unit(c(1, 0.085), "null")
)))
print(normalized_comparison, newpage = FALSE,
vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1))
print(unnormalized_comparison, newpage = FALSE,
vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2))
grid::grid.text(
"Same Y and U: identical accept/reject decisions. Right vertical scale = left / 6.",
gp = grid::gpar(fontsize = 12),
vp = grid::viewport(layout.pos.row = 2, layout.pos.col = 1:2)
)
grid::popViewport()
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.
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, and ?chol. The course’s original reading is Rizzo, M. L. (2007), Statistical Computing with R, Chapman & Hall/CRC.