Review

A review of the whole course: the key algorithms and formulas of Chapters 3–7, and the format of the final exam.

The first two chapters are preliminaries (Chapter 1: discrete-event simulation; Chapter 2: probability and statistics). Chapter 1’s hand simulation of a single-server queue (§1.2) comes back on the exam (Part I, question 3).

Chapter 3: Random Number Generation

Linear congruential generators (LCGs):

$$z_i = (a z_{i-1} + c) \bmod m, \qquad u_i = \frac{z_i}{m} .$$

Necessary and sufficient conditions for a full-period LCG (period $m$, Hull–Dobell theorem). For $c \ne 0$, the LCG has full period if and only if:

1. $c$ and $m$ are relatively prime (their only common divisor is 1);
2. if $q$ is a prime that divides $m$, then $q$ divides $a - 1$;
3. if 4 divides $m$, then 4 divides $a - 1$.

(For a multiplicative LCG, $c = 0$, the period is at most $m - 1$; it is reached when $m$ is prime and $a$ is a primitive element modulo $m$: the PMMLCG.)

Testing RNGs:

1. Test for uniformity: chi-square test. Split $[0, 1]$ into $k$ subintervals of equal length; $f_j$ = number of the $n$ values $u_i$ in interval $j$;

   $$\chi^2 = \frac{k}{n}\sum_{j=1}^k \Big(f_j - \frac{n}{k}\Big)^2 ;$$

   reject uniformity at level $\alpha$ if $\chi^2 > \chi^2_{k-1,\,1-\alpha}$, the $(1-\alpha)$-quantile of $\chi^2_{k-1}$. (The serial test applies the same idea to pairs or $d$-tuples of consecutive $u_i$’s.)
2. Test for independence: runs test (runs up), based on the lengths of the increasing runs in $u_1, u_2, \dots$; the statistic is compared with $\chi^2_6$.

Chapter 4: Generating Random Variates

• Inverse transform method
  1. Generate $U \sim U(0, 1)$.
  2. Return $X = F^{-1}(U)$, where $F$ is the c.d.f.
  • Discrete (or general) case: return $X = \min\{x : F(x) \ge U\}$.
• Composition
  • $F(x) = \sum_{j=1}^\infty p_j F_j(x)$ (a mixture): randomly select one $F_j$ according to the probabilities $p_j$.
  1. Generate a random index $J$ with $P(J = j) = p_j$, $j = 1, 2, \dots$
  2. Return $X$ with distribution $F_J$.
• Convolution: uses the structure of $X$ itself, when $X$ is a sum of i.i.d. r.v.s (e.g. Erlang = sum of exponentials, binomial = sum of Bernoullis).
  • $X = Y_1 + Y_2 + \cdots + Y_m$, where the $Y_i$’s are i.i.d. $\sim G$.
  1. Generate $Y_1, Y_2, \dots, Y_m \sim G$, independent.
  2. Return $X = Y_1 + Y_2 + \cdots + Y_m$.
• Acceptance–rejection: use a surrogate distribution that is easy to sample from to approximate $f$.
  • Majorizing function: $t(x) \ge f(x)$ for all $x$.
  • Normalization: $r(x) = \dfrac{t(x)}{\int_{-\infty}^{\infty} t(x)\, dx} = \dfrac{t(x)}{C}$, where $C = \int t(x)\, dx \ge 1$.
  1. Generate $Y \sim r$.
  2. Generate $U \sim U(0, 1)$, independent of $Y$ ($U \perp Y$).
  3. If $U \le f(Y)/t(Y)$, return $X = Y$; otherwise, go back to step 1 and try again.
  • Efficiency: $P(\text{accept}) = 1/C$, so the expected number of tries per variate is $C$. A good $t$ hugs $f$ closely, making $C$ close to 1.
• Special techniques
  • Box–Muller method (normal)
  • Acceptance–rejection polar algorithm (normal)
  • Table lookup and the alias method: general techniques for discrete distributions (with finitely many values).

Chapter 5: Input Analysis

Goal: select an input distribution to model an input process (e.g. interarrival times, service times). How to do it? Collect data, then three steps:

1. Collect i.i.d. observations of the input random variable of interest (and check independence, §5.1).
2. Step I: perform a statistical analysis (summary statistics, histogram, box plot) and hypothesize a distribution family.
3. Step II: estimate the parameters of the distribution, e.g. by maximum likelihood (MLE):
   1. Construct the likelihood function, the joint density (or p.m.f.) of the data as a function of $\theta$: $L(\theta) = f_\theta(x_1, x_2, \dots, x_n) = \prod_{i=1}^n f_\theta(x_i)$ for i.i.d. data.
   2. MLE: $\hat\theta = \arg\max_\theta L(\theta) = \arg\max_\theta \ln L(\theta)$ (the log turns the product into a sum; set the derivative to 0).
   3. Examples: exponential with mean $\theta$: $\hat\theta = \bar x$; geometric (failures before the first success): $\hat\theta = 1/(1 + \bar x)$; $U(a, b)$: $\hat a = \min x_i$, $\hat b = \max x_i$. An approximate CI for $\theta$ uses the asymptotic normality of the MLE (§5.5).
4. Step III: determine the quality of the fit.
   1. Graphical methods
      1. Q–Q plot (compares quantiles; sensitive in the tails)
      2. P–P plot (compares probabilities; sensitive in the middle; also works for discrete data)
   2. Statistical methods (hypothesis tests of $H_0$: the data are i.i.d. with distribution $\hat F$)
      1. Chi-square test: compare the observed counts $N_j$ in $k$ intervals with the expected counts $np_j$:
         $\chi^2 = \sum_{j=1}^k \frac{(N_j - np_j)^2}{np_j}$; reject if $\chi^2 > \chi^2_{k-1,\,1-\alpha}$.
      2. K–S test: the largest vertical distance between the empirical c.d.f. and $\hat F$:
         $D_n = \max\{D_n^+, D_n^-\}$, $D_n^+ = \max_i\{\frac{i}{n} - \hat F(x_{(i)})\}$, $D_n^- = \max_i\{\hat F(x_{(i)}) - \frac{i-1}{n}\}$; reject if the adjusted statistic exceeds the critical value ($c$, $c'$ or $c''$, by case).
      3. Each test has limitations; given a scenario, you should know which test to apply:

	chi-square test	K–S test
data	discrete or continuous	continuous (the basic form)
sample size	large: every $np_j \ge 5$ (asymptotic test)	any $n$ (exact critical values)
choices to make	the intervals (the result can depend on them; use equal-probability intervals)	none (no grouping of the data)
estimated parameters	degrees of freedom between $k - m - 1$ and $k - 1$	only special cases: all parameters known ($c$), normal ($c'$), exponential ($c''$)

Chapter 6: Output Analysis

Goal: estimate the system (model) performance from simulation outputs.

Idea: construct i.i.d. replications, since the observations within one run are correlated (and not identically distributed).

Terminating simulation (§6.3)

• Run multiple replications. Observations across runs are i.i.d., so average across runs to construct estimators: $X_j$ = output of run $j$, and the CI is

  $$\bar X(n) \pm t_{\alpha/2,\,n-1}\sqrt{\frac{S^2(n)}{n}} .$$

  • Initial conditions: must represent the real system at the start (§6.3.4).
  • Control the estimation precision (§6.3.3):
    • Absolute half-width $\beta$: from $n_0$ pilot runs, the smallest $n$ with $t_{\alpha/2,\,n-1}\sqrt{S^2(n_0)/n} \le \beta$, i.e. roughly $n \approx S^2(n_0)\, t^2_{\alpha/2,\,n-1}/\beta^2$.
    • From half-width $\beta_1$ with $n_1$ runs to $\beta_2$: $n_2 \approx n_1(\beta_1/\beta_2)^2$, i.e. $n_2 - n_1$ more runs.
    • Relative error $\gamma$: the smallest $n$ with $t_{\alpha/2,\,n-1}\sqrt{S^2(n)/n}\,/\,|\bar X(n)| \le \gamma'$, where $\gamma' = \gamma/(1 + \gamma)$ is the adjusted value that guarantees relative error $\gamma$ with respect to the true $\mu$.

Nonterminating simulation (steady-state parameters, §6.4)

• Initial-data deletion (§6.4.1)
  • Delete a number $l$ of observations from the beginning of a run (the warm-up period) and use $\bar Y(m, l) = \frac{1}{m - l}\sum_{i=l+1}^m Y_i$. Trade-off: $l$ too small means bias, $l$ too large means high variance.
  • Use Welch’s procedure to choose $l$: $n$ runs, average across runs, smooth with a moving average of window $w$, plot, and eyeball where the curve levels off.
  • Then replication/deletion (§6.4.2): $n$ runs, each with its first $l$ observations deleted; $X_j = \bar Y_j(m', l)$ are i.i.d., so use the terminating-simulation CI.
• Batch means method (§6.4.3)
  • One run, $Y_1, Y_2, \dots, Y_m$ with $m$ very large.
  • Split it into $n$ batches of size $k$ ($m = nk$), with batch means $\bar Y_1(k), \dots, \bar Y_n(k)$. For large $k$ they are approximately i.i.d. normal, so the CI is $\bar Y(n, k) \pm t_{\alpha/2,\,n-1}\sqrt{S^2_{\bar Y}/n}$. Rule of thumb: $10 \le n \le 30$ batches.
• Regenerative method (§6.4.4)
  • Identify regeneration points, at which the process probabilistically restarts. Observations over different regeneration cycles are i.i.d.: the pairs $(Z_j, N_j)$ = (sum of the observations, number of observations) of cycle $j$.
  • $\mu = \mathbf{E}[Z]/\mathbf{E}[N]$, estimated by $\hat\mu = \bar Z/\bar N$ (biased, but strongly consistent), with CI $\hat\mu \pm z_{\alpha/2}\sqrt{\hat\sigma_V^2/n'}\,/\,\bar N$, $\hat\sigma_V^2 = \hat\sigma_{11} - 2\hat\mu\hat\sigma_{12} + \hat\mu^2\hat\sigma_{22}$.

Multiple performance measures (§6.5)

• $P(V_i \in I_i) = 1 - \alpha_i$, $i = 1, \dots, r$.
• Bonferroni: $P(V_i \in I_i\ \forall i) \ge 1 - \sum_{i=1}^r \alpha_i$, for any dependence between the CIs. E.g. for an overall level of at least 90% with $r = 2$ measures, make two 95% CIs.

Chapter 7: Variance Reduction Techniques

Goal: a smaller $\operatorname{Var}$ for the same number of runs, i.e. a narrower CI at no additional simulation cost.

• Common random numbers (§7.1): to compare systems, $Z = X_1 - X_2$,
  $\operatorname{Var}(Z) = \operatorname{Var}(X_1) + \operatorname{Var}(X_2) - 2\operatorname{cov}(X_1, X_2)$; use the same random numbers for both systems to make $\operatorname{cov} > 0$.
  • Use the inverse transform (monotone, one $u$ per variate).
  • Synchronization: the same $u$ must be used for the same purpose in both systems (dedicated streams, or generate the variates beforehand).
  • Requires both outputs to react to the $u$’s in the same direction; otherwise it backfires.
• Antithetic variates (§7.2): for one system,
  • use the $u_i$’s in one replication and the $1 - u_i$’s in the other;
  • average the pair: $\operatorname{Var}\Big(\frac{X^{(1)} + X^{(2)}}{2}\Big) = \frac{\sigma^2 + \operatorname{cov}(X^{(1)}, X^{(2)})}{2}$, which is better than $\sigma^2/2$ when $\operatorname{cov} < 0$ (output monotone in the $u$’s).
  • Do not combine with CRN.
• Control variates (§7.3)
  • $\mu = \mathbf{E}[X]$; $Y$ can be observed in the same run; $\nu = \mathbf{E}[Y]$ is known; $X$ and $Y$ are correlated.
  • $X_c = X - a(Y - \nu)$ (minus sign), unbiased for every $a$; best $a^* = \operatorname{cov}(X, Y)/\operatorname{Var}(Y)$, giving $\operatorname{Var}(X_c^*) = (1 - \rho_{XY}^2)\operatorname{Var}(X)$.
  • In practice $a^*$ is estimated from the runs, which makes the estimator slightly biased.
• Indirect estimation (§7.4): use known relationships and replace random quantities by known deterministic quantities, e.g. $\tilde w(n) = \hat d(n) + \mathbf{E}[S]$ instead of $\hat w(n)$; with Little’s law ($L = \lambda w$, $L_Q = \lambda d$), also $\tilde L = \lambda\tilde w$ and $\tilde Q = \lambda\hat d$.
• Stratified sampling (§7.5): strata $Y = 1 \mid Y = 2 \mid \cdots$ with known probabilities $P(Y = i)$:
  $\mu = \sum_i \mathbf{E}[X \mid Y = i]\,P(Y = i)$; estimate each $\mathbf{E}[X \mid Y = i]$ by sampling within stratum $i$, and combine with the known weights $P(Y = i)$.
• Conditioning (§7.6): $\mu = \mathbf{E}[X] = \mathbf{E}\big[\mathbf{E}[X \mid Z]\big]$, where $\mathbf{E}[X \mid Z = z]$ can be obtained analytically; use $\mathbf{E}[X \mid Z]$ instead of $X$. Its variance is $\operatorname{Var}(X) - \mathbf{E}[\operatorname{Var}(X \mid Z)] \le \operatorname{Var}(X)$ (law of total variance).
• Importance sampling (§7.7): use a change of measure, so that the desired quantity under the new distribution has a smaller variance:
  $\mu = \mathbf{E}_f[h(Y)] = \mathbf{E}_g\big[h(Y)\,f(Y)/g(Y)\big]$; sample $Y_i \sim g$ and use $\hat\mu = \frac{1}{n}\sum_i h(Y_i)\,f(Y_i)/g(Y_i)$. Choose $g$ roughly proportional to $h f$ (e.g. shifted toward the rare event).

Exam format

Bring your book (for the tables: $z$, $t$, $\chi^2$, K–S critical values) and your calculator.

Part I (3 questions)

1. Random variate generation (Chapter 4): e.g. derive an inverse-transform, composition, convolution or acceptance–rejection algorithm for a given distribution, and compute its efficiency.
2. Input analysis (Chapter 5): MLE (write $L(\theta)$, take the log, differentiate, solve).
3. Queueing simulation (Chapter 1, §1.2). Given a table of interarrival times $A_i$ and service times $S_i$:
   • draw the sample path ($N(t)$ or $Q(t)$, and the server status);
   • calculate performance measures, e.g. the average delay $\hat d(n) = \frac{1}{n}\sum D_i$, the time-average number in queue $\hat q(n) = \frac{1}{T}\int_0^T Q(t)\, dt$, and the server utilization $\hat u(n) = \frac{1}{T}\int_0^T B(t)\, dt$, where $B(t) = 1$ if the server is busy at time $t$ and 0 if idle;
   • fit a distribution to the data (e.g. an exponential, by MLE) and perform a hypothesis test (chi-square or K–S).

Recipe for the hand simulation: arrival times $a_i = a_{i-1} + A_i$; start of service $\max(a_i, d_{i-1})$; delay $D_i = \max(0, d_{i-1} - a_i)$; departure $d_i = a_i + D_i + S_i$. Then list the event times in order, read off $Q(t)$ and $B(t)$ between events, and add up the areas (time length × level).

Part II (3 questions)

1. Output analysis (Chapter 6): construct a CI, and estimate the number of runs needed to achieve a given precision (absolute $\beta$, or relative $\gamma$).
2. Two questions on variance reduction (Chapter 7). For a given performance, measure:
   1. the variance of the crude Monte Carlo estimator (plain sampling, no variance reduction);
   2. the variance of the estimator with the variance reduction technique, and the percentage reduction compared with 1.

example	performance	1. crude variance	2. variance with the technique	reduction
§7.2	$\theta = \mathbf{E}[e^U]$, two evaluations	$0.121$	antithetic: $0.0039$	96.8%
§7.3	$\theta = \mathbf{E}[e^U]$, one evaluation	$0.242$	control variate $U$: $0.0039$	98.4%
§7.5	$\mathbf{E}[e^U]$, $2n$ evaluations	$0.242/2n$	stratified (two halves): $0.065/2n$	73%
§7.7	$P(10 < X < 11)$, $X \sim N(0, 1)$	$\approx 10^{-24}/n$	importance sampling, $g = N(10, 1)$: $\lesssim 10^{-44}/n$	factor $\approx 10^{20}$