Source: Principles of Statistics I, Texas A&M University / Tamhane-Dunlop Ch. 2
Tags: simulation, inverse transform method, quantile function simulation, generating random variables, survival function, hazard function, failure rate, cumulative hazard, reliability, lifetime distribution, R functions, SAS functions
You can simulate observations from any distribution by transforming Uniform(0, 1) draws through the quantile function (the inverse transform method). This works directly when the CDF has a closed-form inverse, and via special algorithms (like Box-Muller for normals) when it does not. In survival and reliability analysis, the survival function S(t), hazard function h(t), and cumulative hazard H(t) provide alternative ways to characterise lifetime distributions, each with a distinct practical interpretation.
Inverse transform method (probability integral transform)
A technique for generating random observations from any distribution. Generate U ~ Uniform(0, 1), then set Y = Q(U) = F⁻¹(U). The resulting Y has CDF F.
Box-Muller algorithm
A method for generating pairs of independent standard normal observations from two Uniform(0, 1) draws. Uses the transformation Z₁ = R cos(2πU₂) and Z₂ = R sin(2πU₂), where R = √(−2 log U₁).
Survival function, S(t)
S(t) = P[T > t] = 1 − F(t). The probability that the event (failure, death) has not yet occurred by time t.
Hazard function (failure rate, intensity function), h(t)
h(t) = f(t)/S(t). The instantaneous rate of failure at time t, given survival up to t. Reported as failures per unit time. For small Δt: Δt · h(t) ≈ P[T ≤ t + Δt | T > t].
Cumulative hazard function, H(t)
H(t) = ∫ from 0 to t of h(τ) dτ. The accumulated instantaneous risk up to time t. Related to survival by S(t) = e^(−H(t)).
Constant failure rate
h(t) = constant for all t. This is the Exponential distribution. Equivalent to the memoryless property.
Increasing failure rate (IFR)
h(t) increases over time. Models wear-out. Weibull with γ > 1 is the standard example.
Decreasing failure rate (DFR)
h(t) decreases over time. Models early-life (infant mortality/burn-in) failures. Weibull with γ < 1.
The core idea is elegant: every continuous CDF F maps its RV to a Uniform(0, 1). Reversing that mapping (using the quantile function) generates observations from F.
The proof
Let Y have strictly increasing continuous CDF F_Y, and let U ~ Uniform(0, 1). Define W = Q_Y(U). Then:
F_W(w) = P[Q_Y(U) ≤ w] = P[U ≤ F_Y(w)] = F_Y(w)
So W has the same distribution as Y.
The recipe
Derive the quantile function Q(u) = F⁻¹(u) in closed form
Generate U from Uniform(0, 1)
Compute Y = Q(U)
Repeat as many times as needed
Worked example: Exponential(β = 4)
CDF: F(y) = 1 − e^(−y/4) for y ≥ 0
Set u = 1 − e^(−yᵤ/4), solve for yᵤ: Q(u) = −4 log(1 − u)
Generate U = runif(1) in R, say U = 0.27
Y = −4 log(1 − 0.27) = −4 log(0.73) ≈ 1.259
Repeat 1000 times for a sample of 1000
Some distributions have CDFs that cannot be inverted in closed form, notably:
Normal distribution (CDF is an integral with no closed-form antiderivative)
Gamma distribution (same issue)
For these, special algorithms exist.
Box-Muller algorithm for Normal(μ, σ²)
Generate U₁, U₂ independently from Uniform(0, 1)
R = √(−2 log U₁)
Z₁ = R cos(2πU₂), Z₂ = R sin(2πU₂)
Z₁ and Z₂ are independent N(0, 1)
Transform: Y₁ = μ + σZ₁, Y₂ = μ + σZ₂ gives independent N(μ, σ²)
This produces two observations per pair of uniform draws.
For a discrete RV with values d₁ < d₂ < … < dₖ and CDF F:
Generate U ~ Uniform(0, 1)
Set D = inf{d : F(d) ≥ U}
In practice: D = dᵢ when F(dᵢ₋₁) < U ≤ F(dᵢ)
Worked example: Geometric(p = 0.2)
PMF: f(i) = 0.2 × 0.8^(i−1), CDF: F(i) = Σ f(k) for k = 1 to i.
F(1) = 0.200
F(2) = 0.360
F(3) = 0.488
F(4) = 0.590
F(5) = 0.672
For U = 0.63: F(4) = 0.590 < 0.63 ≤ 0.672 = F(5), so the simulated value is C = 5.
Proof that this works
f(dⱼ) = P[D = dⱼ] = P[F(dⱼ₋₁) < U ≤ F(dⱼ)] = G(F(dⱼ)) − G(F(dⱼ₋₁)) = F(dⱼ) − F(dⱼ₋₁)
where G is the CDF of U, which is simply G(u) = u on [0, 1].
R uses a consistent naming convention. The prefix letter determines which function you get:
d = density/mass function (PDF or PMF)
p = cumulative distribution function (CDF)
q = quantile function
r = random number generation
Common distributions and their R names (with parameters):
Normal: norm (mean, sd)
Exponential: exp (rate)
Gamma: gamma (shape, rate)
Weibull: weibull (shape, scale)
Beta: beta (shape1, shape2)
Binomial: binom (size, prob)
Poisson: pois (lambda)
Chi-squared: chisq (df)
t: t (df)
F: f (df1, df2)
Uniform: unif (min, max)
Geometric: geom (prob)
Hypergeometric: hyper (m, n, k)
Negative binomial: nbinom (size, prob)
Cauchy: cauchy (location, scale)
Logistic: logis (location, scale)
Lognormal: lnorm (meanlog, sdlog)
Example usage in R:
PDF at y = 4 for N(2, 5): dnorm(4, 2, 5) = 0.0735
CDF at y = 4 for N(2, 5): pnorm(4, 2, 5) = 0.6554
95th percentile for N(2, 5): qnorm(0.95, 2, 5) = 10.224
10 random draws from N(0, 1): rnorm(10, 0, 1)
Critical value for two-sided t-test, 11 df: qt(0.975, 11)
These functions describe the lifetime of a device, organism, or process. Let T be the time to the event of interest, with CDF F and PDF f.
Survival function
S(t) = P[T > t] = 1 − F(t)
The probability the device/subject is still alive (functioning) at time t. Starts at 1, decreases to 0.
Hazard function (failure rate)
h(t) = f(t) / S(t)
The instantaneous risk of failure at time t, given survival up to that point. This is the key function in reliability: it tells you how the risk of failure changes over time.
Constant h(t) = 1/β: Exponential distribution (memoryless)
Increasing h(t): wear-out (e.g. Weibull with γ > 1)
Decreasing h(t): burn-in / infant mortality (e.g. Weibull with γ < 1)
Bathtub-shaped h(t): real devices often have decreasing early failure, constant middle life, increasing late-life failure
Cumulative hazard function
H(t) = ∫ from 0 to t of h(τ) dτ = −log(S(t))
The total accumulated risk up to time t.
Given any one of f(t), S(t), h(t), or H(t), you can derive all the others:
From f(t):
S(t) = ∫ from t to ∞ of f(τ) dτ
h(t) = f(t) / S(t)
H(t) = −log(S(t))
From S(t):
f(t) = −S'(t)
h(t) = −S'(t) / S(t)
H(t) = −log(S(t))
From h(t):
H(t) = ∫ from 0 to t of h(τ) dτ
S(t) = e^(−H(t))
f(t) = h(t) · e^(−H(t))
From H(t):
h(t) = H'(t)
S(t) = e^(−H(t))
f(t) = H'(t) · e^(−H(t))
Exponential(β): h(t) = 1/β (constant), S(t) = e^(−t/β), H(t) = t/β
Weibull(γ, β): h(t) = (γ/β) t^(γ−1) (power of t), S(t) = e^(−t^γ/β)
Gompertz(c, b): h(t) = c · e^(bt) (exponentially increasing hazard). Used in actuarial science and ageing models.
Inverse transform: Y = Q(U) = F⁻¹(U), where U ~ Uniform(0, 1)
Exponential quantile: Q(u) = −β log(1 − u)
Box-Muller: Z₁ = √(−2 log U₁) cos(2πU₂), Z₂ = √(−2 log U₁) sin(2πU₂)
Survival function: S(t) = 1 − F(t) = P[T > t]
Hazard function: h(t) = f(t) / S(t)
Cumulative hazard: H(t) = −log(S(t)) = ∫₀ᵗ h(τ) dτ
Survival from cumulative hazard: S(t) = e^(−H(t))
PDF from hazard: f(t) = h(t) · e^(−∫₀ᵗ h(τ)dτ)
⚠️ The inverse transform method is the foundational simulation technique. Be able to derive Q(u) for a given CDF and walk through the generation steps.
⚠️ Know why the method fails for normal and gamma distributions (no closed-form CDF inverse) and that special algorithms exist.
⚠️ Box-Muller generates two independent normal observations per call. This is efficient and commonly referenced.
⚠️ For discrete simulation, the key step is locating which interval of the CDF the uniform draw falls into. Be comfortable computing cumulative probabilities step by step.
⚠️ Know the four lifetime functions and how to convert between any pair. The conversion table is a likely exam topic.
⚠️ Constant hazard = Exponential = memoryless. If a problem says "constant failure rate," the answer is Exponential.
⚠️ R function naming convention (d/p/q/r prefix) is essential for computational questions. Know what each prefix returns.
Q: Describe the inverse transform method for generating an observation from a distribution with CDF F.
A: Generate U ~ Uniform(0, 1). Compute Y = F⁻¹(U). The observation Y has CDF F. This works when F⁻¹ can be expressed in closed form.
Q: Why can you not use the basic inverse transform method to simulate from a Normal distribution?
A: The normal CDF cannot be inverted in closed form; it is expressed as an integral without an algebraic antiderivative. Special algorithms such as Box-Muller are used instead.
Q: A device has survival function S(t) = e^(−0.05t). What are the hazard function and the distribution of T?
A: H(t) = 0.05t, so h(t) = H'(t) = 0.05 (constant). T has an Exponential distribution with β = 1/0.05 = 20.
Q: Using R, how do you find the probability that a Poisson(25) random variable exceeds 31?
A: 1 − ppois(31, 25), by convention. This gives P[X > 31] = 1 − P[X ≤ 31].
Q: For a Geometric(p = 0.2) distribution, you generate U = 0.45 from a Uniform(0, 1). What is the simulated observation?
A: Compute CDF values: F(1) = 0.2, F(2) = 0.36, F(3) = 0.488. Since F(2) = 0.36 < 0.45 ≤ 0.488 = F(3), the simulated value is 3.
Q: If a Weibull distribution has shape γ = 0.5, what does the hazard function look like and what failure pattern does this represent?
A: The hazard function h(t) is proportional to t^(−0.5), which is decreasing. This models infant mortality or early-life failures, where risk declines as items that were going to fail early are weeded out.
simulation, inverse transform method, probability integral transform, quantile function, random number generation, Box-Muller algorithm, survival function, hazard function, failure rate, intensity function, cumulative hazard, reliability, survival analysis, lifetime distribution, constant failure rate, increasing failure rate, decreasing failure rate, bathtub curve, Exponential hazard, Weibull hazard, Gompertz, R distribution functions, dnorm, pnorm, qnorm, rnorm, SAS RANNOR RANUNI, STAT 301, principles of statistics