Chapter 2 Introduction to Simulation

2.1 Generation of Pseudo-Random Numbers: LCG & Hull-Dobell Theorem

Deterministic computers generate sequences of pseudo-random numbers using mathematical recurrence relations designed to satisfy statistical tests of randomness.

2.1.1 Linear Congruential Generator (LCG)

The LCG generates a sequence of pseudo-random integers \(\{X_n\}\) via: \[X_{n+1} = (a X_n + c) \bmod m\]

\(m\): Modulus (\(m > 0\)).
\(a\): Multiplier (\(0 < a < m\)).
\(c\): Increment (\(0 \le c < m\)).
\(X_0\): Initial seed (\(0 \le X_0 < m\)).

Continuous standard uniform variates \(U_n \in [0, 1)\) are obtained by normalization: \[U_n = \frac{X_n}{m}\]

2.1.1.1 Iterative Calculation Example:

Consider an LCG with parameters \(m = 16, a = 5, c = 3\), and initial seed \(X_0 = 1\):

Iteration 1: \[X_1 = (5 \times 1 + 3) \bmod 16 = 8 \implies U_1 = \frac{8}{16} = 0.5000\]
Iteration 2: \[X_2 = (5 \times 8 + 3) \bmod 16 = 43 \bmod 16 = 11 \implies U_2 = \frac{11}{16} = 0.6875\]
Iteration 3: \[X_3 = (5 \times 11 + 3) \bmod 16 = 58 \bmod 16 = 10 \implies U_3 = \frac{10}{16} = 0.6250\]

2.1.2 Full-Period Conditions: The Hull-Dobell Theorem

An LCG achieves a full period \(m\) (visiting all integers \(0, 1, \dots, m-1\) before repeating) if and only if:

\(c\) and \(m\) are relatively prime (\(\gcd(c, m) = 1\)).
Every prime factor \(p\) dividing \(m\) also divides \(a - 1\) (\(a \equiv 1 \pmod p\)).
If \(4\) divides \(m\), then \(4\) must also divide \(a - 1\) (\(a \equiv 1 \pmod 4\)).

2.2 Generation of Continuous Random Variables

2.2.1 The Inverse-Transform Method

2.2.1.1 Mathematical Formulation:

Let \(U \sim \text{Unif}(0,1)\), and let \(F(x)\) be a strictly increasing continuous Cumulative Distribution Function (CDF). Define \(X = F^{-1}(U)\).

\[\text{Standard Uniform } U \sim \text{Unif}(0,1) \implies \text{Inverse CDF Evaluation } F^{-1}(U) \implies \text{Target Sample } X\]

2.2.1.2 Step-by-Step Mathematical Proof:

We evaluate the CDF of \(X = F^{-1}(U)\):

Step 1 (Definition of CDF):
\[P(X \le x) = P(F^{-1}(U) \le x)\]
Step 2 (Monotonic Transformation):
Rationale: Because \(F(x)\) is a strictly increasing function, applying \(F\) to both sides of an inequality preserves the inequality direction:
\[P(F(F^{-1}(U)) \le F(x)) = P(U \le F(x))\]
Step 3 (Uniform Distribution Property):
Rationale: For any standard uniform random variable \(U \sim \text{Unif}(0,1)\), \(P(U \le u) = u\) for all \(u \in [0, 1]\). Setting \(u = F(x)\):
\[P(U \le F(x)) = F(x)\]
Conclusion: \(P(X \le x) = F(x)\), proving that \(X = F^{-1}(U)\) follows CDF \(F(x)\). \(\blacksquare\)

2.2.1.3 Application: Generating Exponential Variates (\(X \sim \text{Exp}(\lambda)\))

Target CDF: \(F(x) = 1 - e^{-\lambda x}\) for \(x \ge 0\).
Set \(U = 1 - e^{-\lambda X}\).
Solve for \(X\):
\[1 - U = e^{-\lambda X} \implies \ln(1 - U) = -\lambda X \implies X = -\frac{\ln(1 - U)}{\lambda}\]
Since \(1 - U \sim \text{Unif}(0,1)\), this simplifies to: \[\mathbf{X = -\frac{\ln(U)}{\lambda}}\]

2.2.2 The Accept-Reject Algorithm

2.2.2.1 Mathematical Formulation:

When \(F^{-1}(u)\) lacks a closed-form expression, we sample candidate variates \(Y\) from a proposal density \(g(y)\) that bounds target density \(f(y)\) under an envelope \(c \cdot g(y)\) (\(c \ge 1\), \(\forall y: f(y) \le c \cdot g(y)\)).

\[\text{Target Density } f(y) \le \text{Envelope Density } c \cdot g(y) \quad \forall y\]

2.2.2.2 Algorithm Steps:

Draw candidate \(Y \sim g(y)\) and independent uniform variate \(U \sim \text{Unif}(0,1)\).
Accept \(X = Y\) if \(U \le \frac{f(Y)}{c \cdot g(Y)}\); otherwise reject and repeat.

2.2.2.3 Proof of Correctness:

Step 1 (Joint Probability of Candidate & Acceptance):
\[P(Y \le x, \text{Accept}) = \int_{-\infty}^x \left[ \int_0^{\frac{f(y)}{c g(y)}} 1 \, du \right] g(y) \, dy = \int_{-\infty}^x \frac{f(y)}{c g(y)} g(y) \, dy = \frac{1}{c} \int_{-\infty}^x f(y) \, dy = \frac{F(x)}{c}\]
Step 2 (Overall Acceptance Probability):
Rationale: Setting \(x = \infty\) evaluates the marginal probability of acceptance per iteration:
\[P(\text{Accept}) = \frac{1}{c} \int_{-\infty}^\infty f(y) \, dy = \frac{1}{c}\]
Step 3 (Conditional Distribution):
Rationale: Applying Bayes’ rule for conditional probability \(P(A \mid B) = \frac{P(A \cap B)}{P(B)}\):
\[P(Y \le x \mid \text{Accept}) = \frac{P(Y \le x, \text{Accept})}{P(\text{Accept})} = \frac{F(x)/c}{1/c} = F(x)\]
Conclusion: Accepted variates follow target density \(f(x)\) with an expected number of draws per accepted sample equal to \(c\). \(\blacksquare\)

2.2.3 Generating Normal Vectors: Cholesky Decomposition

To generate a correlated multivariate Gaussian vector \(X \sim \mathcal{N}_n(\mu, \Sigma)\):

Compute the Cholesky Factorization of covariance matrix \(\Sigma\):
\[\Sigma = L L^T \quad \text{(where } L \text{ is lower-triangular)}\]
Draw independent standard normal vector \(Z = (Z_1, \dots, Z_n)^T \sim \mathcal{N}_n(0, I_n)\).
Compute \(X = \mu + L Z\).

2.2.3.1 Proof of Covariance Structure:

\[\text{Cov}(X) = \text{Cov}(\mu + L Z) = L \text{Cov}(Z) L^T = L I_n L^T = L L^T = \Sigma \quad \blacksquare\]

2.3 Simulation of Stochastic Processes

2.3.1 Homogeneous Poisson Process (HPP)

In an HPP with constant rate \(\lambda\), inter-arrival times between consecutive events are independent exponential variates \(T_k \sim \text{Exp}(\lambda)\).

Algorithm:
Set \(S_0 = 0\). For \(k = 1, 2, \dots\): draw \(U_k \sim \text{Unif}(0,1)\), compute \(T_k = -\frac{\ln U_k}{\lambda}\), and set event arrival time \(S_k = S_{k-1} + T_k\).

2.3.2 Non-Homogeneous Poisson Process (NHPP): Thinning Method

When arrival intensity \(\lambda(t)\) varies over time but is bounded by \(\lambda_{\max}\):

Generate candidate arrival times \(S_k\) from an HPP with constant rate \(\lambda_{\max}\).
Accept candidate \(S_k\) with probability \(\frac{\lambda(S_k)}{\lambda_{\max}}\) (draw \(U \sim \text{Unif}(0,1)\), accept if \(U \le \frac{\lambda(S_k)}{\lambda_{\max}}\)).

2.4 Random Permutations: Fisher-Yates Shuffle

To generate an unbiased random permutation of an array of \(n\) elements such that all \(n!\) permutations are equiprobable:

Algorithm:
For \(i = n, n-1, \dots, 2\):
1. Draw a random integer \(j\) uniformly from \(\{1, 2, \dots, i\}\).
2. Swap elements at index \(i\) and index \(j\).

2.5 Empirical Distribution Function (EDF) & Influence Functions

2.5.1 Empirical Distribution Function (EDF)

Given an observed sample \(X_1, \dots, X_n\), the EDF \(F_n(x)\) represents the empirical step CDF: \[F_n(x) = \frac{1}{n} \sum_{i=1}^n \mathbb{I}(X_i \le x)\]

2.5.1.1 Glivenko-Cantelli Theorem:

As sample size \(n \to \infty\), \(F_n(x)\) converges uniformly to the true underlying CDF \(F(x)\) almost surely: \[\sup_{x \in \mathbb{R}} |F_n(x) - F(x)| \xrightarrow{a.s.} 0\]

2.5.2 Influence Function (IF)

The Influence Function quantifies the sensitivity of a functional estimator \(T(F)\) to an infinitesimal point-mass contamination at \(x\): \[IF(x; T, F) = \lim_{\varepsilon \to 0} \frac{T((1-\varepsilon)F + \varepsilon \delta_x) - T(F)}{\varepsilon}\] For the sample mean, \(IF(x; \text{Mean}, F) = x - \mu\) (demonstrating linear, unbounded sensitivity to extreme outliers).