7. Variance Reduction Techniques

Ways to get tighter confidence intervals from the same simulation budget: common random numbers, antithetic variates, control variates, indirect estimation, stratified sampling, conditioning, and importance sampling.

Estimate $\mu = \mathbf{E}[X]$ from $n$ i.i.d. replications $X_1, \dots, X_n$: the point estimator is $\hat\mu = \bar X(n)$, and the $100(1-\alpha)\%$ CI is

$$\left[\hat\mu - t_{\alpha/2,\,n-1}\sqrt{\frac{\sigma^2}{n}},\;\; \hat\mu + t_{\alpha/2,\,n-1}\sqrt{\frac{\sigma^2}{n}}\,\right], \qquad \sigma^2 = \operatorname{Var}(X)$$

(in practice $\sigma^2$ is estimated by $S^2(n)$, as in §6.3).

Goal: reduce the variance $\sigma^2$ at no additional simulation cost. The half-width is proportional to $\sigma/\sqrt n$, so there are two ways to make the CI narrower: more replications (a larger $n$, which costs computer time; by §6.3.3 the required $n$ is proportional to $\sigma^2$), or a smaller $\sigma^2$ for the same $n$. Variance reduction techniques do the second, by using the random numbers more cleverly: halving $\sigma^2$ is as good as doubling $n$, but costs (almost) nothing.

7.1. Common Random Numbers (CRN)

Applicable when comparing two or more alternative systems.

Idea: we want the simulation experiments to be conducted under the same “experimental conditions”, i.e. with the same random numbers, so that an observed difference in performance is due to the difference between the systems and not to the random numbers (“luck”) of each simulation.

Ex: 1. M/M/1, 2. M/M/2.

$X_i$ = average delay in queue in system $i = 1, 2$. Goal: estimate $\zeta = \mathbf{E}[Z] = \mathbf{E}[X_1 - X_2]$, where $Z = X_1 - X_2$.

Repeat the experiment $n$ times (each time simulating both systems) and get $Z_1, Z_2, \dots, Z_n$, with $Z_j = X_{1j} - X_{2j}$. The point estimator is

$$\bar Z(n) = \frac{1}{n}\sum_{j=1}^n Z_j ,$$

and the half-width of the $100(1-\alpha)\%$ CI is

$$t_{\alpha/2,\,n-1}\sqrt{\frac{\operatorname{Var}(Z_1)}{n}}$$

(with $\operatorname{Var}(Z_1)$ estimated by the sample variance of the $Z_j$’s), where

$$\operatorname{Var}(Z_1) = \operatorname{Var}(X_1 - X_2) = \operatorname{Var}(X_1) + \operatorname{Var}(X_2) - 2\operatorname{cov}(X_1, X_2) .$$

If the two systems are simulated with independent random numbers, $\operatorname{cov}(X_1, X_2) = 0$. We want $\operatorname{cov}(X_1, X_2) > 0$: then the last term is subtracted and $\operatorname{Var}(Z_1)$ is smaller than with independent sampling.

CRN is a technique that tries to induce this positive correlation by using the same random numbers to simulate all alternative systems. (A first example was in §6.3.3: for M/M/1 vs M/D/1, CRN reduced $\operatorname{Var}(Z_j)$ from 3.79 to 2.23.)

When does CRN work?

$$u_1, u_2, \dots, u_n \;\longrightarrow\; \boxed{\text{simulation program}} \;\longrightarrow\; X$$

The simulation program is a (complicated, but deterministic) function of the random numbers: $X = X(u_1, u_2, \dots)$. To picture it, pretend that $X$ depends on a single random number: $X = X(u)$. CRN feeds the same $u$ to both systems, so the question is how $X_1(u)$ and $X_2(u)$ react to $u$.

Desired cases: $X_1$ and $X_2$ both increase, or both decrease, in $u$.

Backfire: one increases while the other decreases. Then CRN induces a negative correlation and increases the variance.

How can we guarantee the desired case? There is no guarantee: the simulation program is far too complicated to check whether it is monotone in every random number. But we can make it very likely, by looking at the two parts that turn random numbers into an output:

flowchart LR
    U(["uᵢ's"]) --> RVG
    subgraph H["X(u₁, u₂, …): monotone if each box is monotone"]
        RVG["RVG: random numbers → random variates"] --> SP["simulation program: random variates → X"]
    end
    SP --> X(["X"])

The composition of two monotone functions is also monotone, so each part has to be monotone:

• The random variate generator (RVG) is under our control. Use the inverse transform (§4.2.1): $X = F^{-1}(u)$ is increasing in $u$, so a larger $u$ always gives a larger variate. Other methods are not monotone; e.g. with acceptance–rejection a slightly different $u$ can change which candidate is accepted, and even how many $u$’s are used. So always use the inverse transform in the CRN approach.
• The simulation program (the model itself) should react monotonically to the variates. This we cannot control, but it can usually be checked by reasoning about the model, as in the queueing example below.

Queueing simulation (service times in the M/M/1 and the M/M/2 system, both generated by the inverse transform from the same $u$):

$$u \uparrow \;\Rightarrow\; S_1 \uparrow \text{ and } S_2 \uparrow \;\Rightarrow\; X_1 \uparrow \text{ and } X_2 \uparrow$$

(longer service times mean longer delays in both systems). Likewise for interarrival times: $u \uparrow \Rightarrow A_1 \uparrow$ and $A_2 \uparrow \Rightarrow X_1 \downarrow$ and $X_2 \downarrow$. Both systems react in the same direction to every random number, which is the desired case.

Synchronization: it is not enough to use the same stream of random numbers to simulate the alternative systems. We must match up (synchronize) the random numbers across the different systems, to ensure that the same $u$ is used for the same purpose in the different systems.

Ex: one stream $u_1, u_2, \dots$ is used for everything, each random number being taken when the program needs one (A = used for an interarrival time, S = used for a service time):

random number	used in M/M/1 for	used in M/M/2 for
$u_1$	A	A
$u_2$	A	A
$\vdots$	$\vdots$	$\vdots$
$u_7$	S	S
$u_8$	A	A
$\to u_9$	A	S

Up to $u_8$ the two systems happen to need the random numbers for the same purposes. But the systems are different, so their events eventually happen in a different order: e.g. in the M/M/2 system a waiting customer starts service on the second server (which needs a service time) while in the M/M/1 system the next event is an arrival (which needs an interarrival time). From $u_9$ on, the same $u$ is used as an interarrival time in one system and as a service time in the other, and all later random numbers are shifted as well. The two systems are no longer driven by the same inputs, and most of the positive correlation is lost.

General synchronization methods:

• Dedicate a stream of random numbers to a particular source of randomness (§3.3), e.g. stream 1 for all interarrival times and stream 2 for all service times. Then the $k$-th interarrival time always uses the $k$-th number of stream 1, in both systems, no matter in which order the events happen.
• Use the inverse transform method to generate random variates: it guarantees monotonicity and is one-to-one (exactly one $u$ per variate). With acceptance–rejection, the number of $u$’s used per variate is random, so the streams drift out of synchronization even when they are dedicated.
• Generate all random variates beforehand, e.g. generate customer $j$’s service time when customer $j$ arrives and attach it to the customer. Then customer $j$ has the same service time in both systems, even if it starts service at different times or (with other queue disciplines) in a different order.

E.g. with numbers (my choice): arrival rate 1; system 1 has one server with mean service time 0.9, system 2 has two servers with mean service time 1.8 each (the same total capacity); $X_i$ = average delay in queue of the first 100 customers, starting empty; 3000 replications of each version:

random numbers	$\operatorname{corr}(X_1, X_2)$	$\operatorname{Var}(Z)$	change
independent for the two systems	0.01	16.4
CRN, one shared stream, not synchronized	0.08	15.7	−4%
CRN, dedicated streams for arrivals and services (synchronized)	0.997	0.057	−99.7%
backfire: system 2 uses $1 - u$ wherever system 1 uses $u$	−0.36	23.3	+42%

Here $\mathbf{E}[Z] \approx 0.43$ (the M/M/2 system has slightly shorter delays). With $n = 10$ replications, the 95% CI half-width is $2.262\sqrt{16.4/10} = 2.9$ with independent random numbers (the CI cannot even tell which system is better), but $2.262\sqrt{0.057/10} = 0.17$ with synchronized CRN. Note that sharing a stream without synchronization gives almost nothing: CRN only works when the same $u$ is used for the same purpose.

7.2. Antithetic Variates (AV)

Used when investigating the performance of a single system.

Idea: make the replications in pairs, such that a large observation in the first replication of a pair tends to be matched by a small one in the second; their average is then closer to $\mu$ than the average of two independent replications.

• Use $u_i$’s for the first replication of a pair to get $X^{(1)}$.
• Use $1 - u_i$’s for the second replication, in the same roles, to get $X^{(2)}$.

Since $1 - U \sim U(0, 1)$ when $U \sim U(0, 1)$, both are valid replications: $\mathbf{E}[X^{(1)}] = \mathbf{E}[X^{(2)}] = \mu$ and $\operatorname{Var}(X^{(1)}) = \operatorname{Var}(X^{(2)}) = \sigma^2$. ($X^{(1)}, X^{(2)}$ are the two runs of one pair, not two systems as in §7.1.) Use the pair average as the observation:

$$X = \frac{X^{(1)} + X^{(2)}}{2}, \qquad \mathbf{E}[X] = \mu,$$

$$\operatorname{Var}(X) = \frac{1}{4}\Big[\operatorname{Var}(X^{(1)}) + \operatorname{Var}(X^{(2)}) + 2\operatorname{cov}(X^{(1)}, X^{(2)})\Big] = \frac{\sigma^2 + \operatorname{cov}(X^{(1)}, X^{(2)})}{2} .$$

With two independent replications the covariance would be 0 and the variance $\sigma^2/2$. So AV helps exactly when $\operatorname{cov}(X^{(1)}, X^{(2)}) < 0$. With $n$ pairs (different pairs use independent random numbers), the pair averages $X_1, \dots, X_n$ are i.i.d., and the CI is computed from them as usual, with $t_{\alpha/2,\,n-1}$.

Why the covariance is negative:

$$u_i \uparrow \;\Rightarrow\; X^{(1)} \uparrow, \qquad 1 - u_i \downarrow \;\Rightarrow\; X^{(2)} \downarrow .$$

Applicability: the output reacts monotonically to the underlying random numbers (the same condition as for CRN; here it makes the two runs of a pair move in opposite directions).

Practical issues:

• Synchronization: $1 - u_i$ must be used for exactly the same purpose in the second run as $u_i$ in the first.
• Use the inverse transform to generate random variates (monotone, one $u$ per variate).
• Use dedicated streams for the different sources of randomness.
• Do not combine CRN with AV (in a comparison of systems): it may backfire. With both, system 1’s first run (with $u$) and system 2’s second run (with $1 - u$) are negatively correlated, and these cross-covariances enter $\operatorname{Var}(Z)$ with the wrong sign: they increase it. Whether the net effect is a reduction is not clear in advance.

Ex: $\theta = \mathbf{E}[e^U]$, $U \sim U(0, 1)$:

$$\theta = \int_0^1 e^u \cdot 1\, du = e - 1 = 1.718 .$$

To compare fairly, both approaches use two evaluations of $e^u$.

Approach 1: crude Monte Carlo simulation. Generate $U_1, U_2 \sim U(0, 1)$, independent ($U_1 \perp U_2$), and use

$$\bar\theta = \frac{1}{2}\left(e^{U_1} + e^{U_2}\right).$$

$$\operatorname{Var}(\bar\theta) = \frac{1}{4}\cdot 2\operatorname{Var}(e^{U}) = \frac{\operatorname{Var}(e^U)}{2} = \frac{\mathbf{E}[e^{2U}] - (e - 1)^2}{2},$$

where

$$\mathbf{E}[e^{2U}] = \int_0^1 e^{2u}\, du = \frac{e^2 - 1}{2},$$

so

$$\operatorname{Var}(\bar\theta) = \frac{1}{2}\left[\frac{e^2 - 1}{2} - (e - 1)^2\right] = \frac{0.2420}{2} = 0.121 .$$

Approach 2: AV. Generate one $U$ and use

$$\hat\theta = \frac{1}{2}\left(e^U + e^{1-U}\right).$$

$$\begin{aligned} \operatorname{Var}(\hat\theta) &= \frac{1}{4}\Big[\operatorname{Var}(e^U) + \operatorname{Var}(e^{1-U}) + 2\operatorname{cov}(e^U, e^{1-U})\Big] \\ &= \frac{1}{2}\operatorname{Var}(e^U) + \frac{1}{2}\operatorname{cov}(e^U, e^{1-U}) = 0.121 + \frac{1}{2}\operatorname{cov}(e^U, e^{1-U}), \end{aligned}$$

using $\operatorname{Var}(e^{1-U}) = \operatorname{Var}(e^U)$ (since $1 - U \sim U(0,1)$). The covariance:

$$\operatorname{cov}(e^U, e^{1-U}) = \mathbf{E}\big[e^U e^{1-U}\big] - \mathbf{E}[e^U]\,\mathbf{E}[e^{1-U}] = e - (e - 1)^2 = -0.2342,$$

because $e^U e^{1-U} = e$ is a constant and $\mathbf{E}[e^{1-U}] = \mathbf{E}[e^U] = e - 1$. So

$$\operatorname{Var}(\hat\theta) = 0.121 - \frac{0.2342}{2} = 0.0039 .$$

The percentage reduction in variance is $(0.121 - 0.0039)/0.121 = 96.8\%$.

Why so much: $e^u$ is increasing in $u$, so $e^U$ and $e^{1-U}$ move in opposite directions, and their average is almost constant:

In a real simulation the reduction is usually much smaller, because the output is far from such a simple, almost linear function of each $u$. E.g. the M/M/1 queue above ($X$ = average delay of the first 100 customers, $\mu \approx 4.13$, $\sigma^2 \approx 8.67$), with $1 - u$ for both the interarrival and the service times in the second run: $\operatorname{corr}(X^{(1)}, X^{(2)}) = -0.36$, so the pair average has variance $2.80$ instead of $\sigma^2/2 = 4.34$, a reduction of 35%.

7.3. Control Variates (CV)

Goal: estimate $\mu = \mathbf{E}[X]$, $X$ = average delay in queue of the first 100 customers.

Suppose there is an intermediate r.v. $Y$ satisfying

1. the value of $Y$ can be observed (in the same run as $X$);
2. $Y$ and $X$ are correlated;
3. $\nu = \mathbf{E}[Y]$ is known.

E.g. $Y$ = average service time of the first 99 customers. $\nu$ is known: it is the mean of the service-time distribution, which we chose as an input. (99, not 100: the delays of customers $1, \dots, 100$ depend only on the service times of customers $1, \dots, 99$.)

If $Y > \nu$ (the service times happened to be long in this run) $\Rightarrow$ probably $X > \mu$ (the delays are long too) $\Rightarrow$ adjust the value of $X$ downwards, towards $\mu$.

If $Y < \nu$ $\Rightarrow$ probably $X < \mu$ $\Rightarrow$ adjust $X$ upwards.

$Y$ is called a control variate: it measures how lucky or unlucky the run was, and we correct $X$ for that luck.

How to decide the adjustment amount?

Let $X$ = the observed value and $X_c$ = the value of $X$ after the adjustment, so the adjustment amount is $X - X_c$. Assume it is proportional to how far $Y$ is from its mean:

$$X - X_c = a(Y - \nu), \qquad X_c = X - a(Y - \nu).$$

$$\mathbf{E}[X_c] = \mathbf{E}[X] - a\,\mathbf{E}[Y - \nu] = \mu - a \cdot 0 = \mu,$$

so $X_c$ is unbiased for every constant $a$. Its variance ($\nu$ is a constant, so it does not affect the variance):

$$\operatorname{Var}(X_c) = \operatorname{Var}(X - aY) = \operatorname{Var}(X) + \underbrace{\big\{a^2\operatorname{Var}(Y) - 2a\operatorname{cov}(X, Y)\big\}}_{\text{want } \le\, 0} .$$

$\operatorname{Var}(X_c) \le \operatorname{Var}(X)$ exactly when the term in braces is $\le 0$.

If $a = 1$, then $Y$ should be chosen such that $\operatorname{Var}(Y) \le 2\operatorname{cov}(X, Y)$.

Idea: find the $a^*$ such that $\operatorname{Var}(X_c)$ is minimized:

$$\frac{d\operatorname{Var}(X_c)}{da} = 2a\operatorname{Var}(Y) - 2\operatorname{cov}(X, Y) = 0 \quad\Rightarrow\quad a^* = \frac{\operatorname{cov}(X, Y)}{\operatorname{Var}(Y)}$$

(a minimum, since the second derivative $2\operatorname{Var}(Y) > 0$). Plugging $a^*$ in:

$$\operatorname{Var}(X_c^*) = \operatorname{Var}(X) - \frac{\operatorname{cov}(X, Y)^2}{\operatorname{Var}(Y)} = \big(1 - \rho_{XY}^2\big)\operatorname{Var}(X), \qquad \rho_{XY} = \operatorname{corr}(X, Y).$$

So with the best $a$, the variance is reduced by the fraction $\rho_{XY}^2$: the more strongly $X$ and $Y$ are correlated (positively or negatively), the better. $a^*$ has the sign of $\operatorname{cov}(X, Y)$: with a negatively correlated control, $a^* < 0$ and $X$ is adjusted upwards when $Y > \nu$.

E.g. M/M/1 queue (arrival rate 1, mean service time 0.9, so $\nu = 0.9$), $X$ and $Y$ as above: $\operatorname{Var}(X) = 8.67$, $\operatorname{Var}(Y) = 0.00819$, $\operatorname{cov}(X, Y) = 0.1663$, so $\rho_{XY} = 0.624$ and $a^* = 0.1663/0.00819 = 20.3$. Then $\operatorname{Var}(X_c^*) = (1 - 0.624^2) \times 8.67 = 5.29$, a reduction of 39%. With $a = 1$ instead, the condition $\operatorname{Var}(Y) \le 2\operatorname{cov}(X, Y)$ holds ($0.008 \le 0.33$), but the variance only drops to 8.35 (4%): $Y$ varies much less than $X$, so a large factor $a$ is needed.

One limitation: $\operatorname{Var}(Y)$ and $\operatorname{cov}(X, Y)$, and hence $a^*$, are usually unknown. In practice, perform $n$ simulation runs, giving $(X_1, Y_1), \dots, (X_n, Y_n)$, and estimate them from the runs:

$$s_Y^2(n) = \frac{\sum_{i=1}^n \big(Y_i - \bar Y(n)\big)^2}{n-1} \approx \operatorname{Var}(Y),$$

$$\hat C_{XY}(n) = \frac{\sum_{i=1}^n \big(X_i - \bar X(n)\big)\big(Y_i - \bar Y(n)\big)}{n-1} \approx \operatorname{cov}(X, Y),$$

$$\hat a^*(n) = \frac{\hat C_{XY}(n)}{s_Y^2(n)}$$

(the slope of the least-squares line of $X$ on $Y$, the line in the figure above). The control-variate estimator is then

$$\bar X_c(n) = \bar X(n) - \hat a^*(n)\big[\bar Y(n) - \nu\big].$$

Drawback: $\hat a^*(n)$ is now a random variable computed from the same $Y_i$’s as $\bar Y(n)$, so the two are not independent, and

$$\mathbf{E}\Big[\hat a^*(n)\big(\bar Y(n) - \nu\big)\Big] \;\ne\; \mathbf{E}\big[\hat a^*(n)\big]\,\mathbf{E}\big[\bar Y(n) - \nu\big] = 0 .$$

So $\mathbf{E}[\bar X_c(n)] \ne \mu$: the estimator is biased, and the usual CI formula is no longer exact. (This is the standard drawback; the formula on the board was most likely this one.) The bias goes to 0 as $n$ grows. It can be avoided by estimating $a^*$ from separate pilot runs: then $\hat a^*$ is independent of $\bar Y(n)$, the expectation above factorizes, and it equals 0. E.g. in the M/M/1 example with $n = 10$ runs and $\hat a^*$ estimated from them, the variance of $\bar X_c(10)$ was 0.58 instead of 0.86 for $\bar X(10)$ (a 33% reduction), with a bias of $-0.07$.

Ex. Estimate $\theta = \mathbf{E}[e^U]$, where $U \sim U(0, 1)$ ($\theta = e - 1 = 1.718$, as in §7.2).

Crude Monte Carlo: generate $U \sim U(0, 1)$ and use $\tilde\theta = e^U$, with

$$\operatorname{Var}(\tilde\theta) = \operatorname{Var}(e^U) = 0.242 .$$

Control variates: let $U$ itself be the control variate. It is observed, it is correlated with $e^U$, and its mean $\mathbf{E}[U] = \frac{1}{2}$ is known.

$$\hat\theta = e^U - a^*\left(U - \tfrac{1}{2}\right), \qquad a^* = \frac{\operatorname{cov}(e^U, U)}{\operatorname{Var}(U)} = 1.690,$$

because

$$\operatorname{cov}(e^U, U) = \mathbf{E}[Ue^U] - \mathbf{E}[U]\,\mathbf{E}[e^U] = 1 - \tfrac{1}{2}(e - 1) = 0.1409,$$

where

$$\mathbf{E}[Ue^U] = \int_0^1 u e^u\, du = \Big[(u - 1)e^u\Big]_0^1 = 1,$$

and $\operatorname{Var}(U) = \frac{1}{12}$, so $a^* = 12 \times 0.1409 = 1.690$ (1.692 on the board is a small rounding slip). So

$$\hat\theta = e^U - 1.690\left(U - \tfrac{1}{2}\right),$$

$$\operatorname{Var}(\hat\theta) = \operatorname{Var}(\tilde\theta)\left(1 - \frac{\operatorname{cov}^2(e^U, U)}{\operatorname{Var}(e^U)\operatorname{Var}(U)}\right) = 0.242 \times (1 - 0.9837) = 0.0039$$

(this is $(1 - \rho^2)\operatorname{Var}(X)$ from above, with $\rho^2 = 0.1409^2/(0.242 \times \frac{1}{12}) = 0.9837$).

Reduction in variance: $(0.242 - 0.0039)/0.242 = 98.4\%$.

Why so much: $e^u$ is almost a straight line on $[0, 1]$, and $1.718 + 1.690(u - \frac{1}{2})$ is the best straight line through it, so $\hat\theta$ only keeps the small curvature. (Compare §7.2: antithetic variates also reached 0.0039, but with two evaluations of $e^u$; here one is enough.)

7.4. Indirect Estimation

Idea: use known relationships, and replace random quantities by constant known quantities wherever possible.

Queueing example:

• $D_i$ = delay in queue of the $i$-th customer
• $W_i$ = total waiting time in system of the $i$-th customer
• $S_i$ = service time of the $i$-th customer
• $Q(t)$ = number of customers in queue at time $t$
• $L(t)$ = number of customers in system at time $t$

Quantities of interest (steady state):

$$d = \lim_{i \to \infty}\mathbf{E}[D_i] = \lim_{n \to \infty}\frac{1}{n}\sum_{i=1}^n D_i, \qquad w = \lim_{i \to \infty}\mathbf{E}[W_i] = \lim_{n \to \infty}\frac{1}{n}\sum_{i=1}^n W_i,$$

$$L = \lim_{t \to \infty}\mathbf{E}[L(t)], \qquad L_Q = \lim_{t \to \infty}\mathbf{E}[Q(t)]$$

(the second equalities hold with probability 1, by the same long-run-average argument as in Chapter 6, §6.4.4). Their direct estimators, from a run with $n$ customers that ends at time $T(n)$ (e.g. when the $n$-th customer leaves):

$$\hat d(n) = \frac{1}{n}\sum_{i=1}^n D_i, \qquad \hat w(n) = \frac{1}{n}\sum_{i=1}^n W_i,$$

$$\hat L(n) = \frac{1}{T(n)}\int_0^{T(n)} L(t)\, dt, \qquad \hat Q(n) = \frac{1}{T(n)}\int_0^{T(n)} Q(t)\, dt$$

(the last two are time averages, as in Chapter 1).

Since $W_i = D_i + S_i$,

$$\hat w(n) = \frac{1}{n}\sum_{i=1}^n (D_i + S_i) = \frac{1}{n}\sum_{i=1}^n D_i + \underbrace{\left\{\frac{1}{n}\sum_{i=1}^n S_i\right\}}_{\to\ \mathbf{E}[S]} .$$

The term in braces is a random average of the service times, but its mean $\mathbf{E}[S]$ is known (we chose the service-time distribution). Replace it by its known mean:

$$\tilde w(n) = \frac{1}{n}\sum_{i=1}^n D_i + \mathbf{E}[S] = \hat d(n) + \mathbf{E}[S] .$$

Both have the same expectation. And

$$\operatorname{Var}(\hat w(n)) = \operatorname{Var}(\hat d(n)) + \operatorname{Var}(\bar S(n)) + 2\operatorname{cov}(\hat d(n), \bar S(n)), \qquad \operatorname{Var}(\tilde w(n)) = \operatorname{Var}(\hat d(n)),$$

where $\bar S(n) = \frac{1}{n}\sum_i S_i$. In a queue, long service times cause long delays, so $\operatorname{cov}(\hat d(n), \bar S(n)) > 0$ and $\operatorname{Var}(\tilde w(n)) < \operatorname{Var}(\hat w(n))$. ($\tilde w(n) = \hat w(n) - [\bar S(n) - \mathbf{E}[S]]$ is a control-variate estimator with $Y = \bar S(n)$ and $a = 1$, §7.3.)

Little’s Law

$$L = \lambda \cdot w$$

where $\lambda$ is the arrival rate; likewise $L_Q = \lambda \cdot d$.

Why: each customer is charged 1 dollar for every unit of time he spends in the system. Compute the long-run average reward earned per unit time in two ways:

• By time: at time $t$, $L(t)$ customers are in the system, each paying 1 dollar per unit time, so the reward rate is $L(t)$. The total reward up to time $T$ is $\int_0^T L(t)\, dt$, and the average reward per unit time is $\frac{1}{T}\int_0^T L(t)\, dt \to L$.
• By customer: customer $i$ pays $W_i$ in total. With $N(T)$ = the number of customers who arrived by time $T$,

$$\frac{1}{T}\sum_{i=1}^{N(T)} W_i = \frac{N(T)}{T}\cdot\frac{1}{N(T)}\sum_{i=1}^{N(T)} W_i \;\to\; \lambda \cdot w,$$

since $N(T)/T \to \lambda$ (arrivals per unit time) and the average of the $W_i$’s $\to w$. (Customers still in the system at time $T$ have not paid everything yet, but this is negligible as $T \to \infty$.)

Both compute the same long-run reward rate, so $L = \lambda w$. Charging only for the time spent in the queue gives $L_Q = \lambda d$ in the same way.

Use in indirect estimation: $\lambda$ is a known input, so $L$ and $L_Q$ can also be estimated indirectly, instead of by the time averages $\hat L(n)$, $\hat Q(n)$:

$$\tilde L(n) = \lambda\,\tilde w(n) = \lambda\big(\hat d(n) + \mathbf{E}[S]\big), \qquad \tilde Q(n) = \lambda\,\hat d(n) .$$

So all four quantities are estimated from $\hat d(n)$ alone.

E.g. M/M/1 queue ($\lambda = 1$, $\mathbf{E}[S] = 0.9$, empty start), first $n = 1000$ customers, 20,000 runs:

quantity	direct estimator	variance	indirect estimator	variance	reduction
$w$	$\hat w(n)$	14.48	$\tilde w(n) = \hat d(n) + \mathbf{E}[S]$	14.37	0.8%
$L$	$\hat L(n)$	16.05	$\tilde L(n) = \lambda\tilde w(n)$	14.37	10.5%
$L_Q$	$\hat Q(n)$	15.83	$\tilde Q(n) = \lambda\hat d(n)$	14.37	9.3%

The gains are small here (the delays, $d \approx 8.1$, dominate $W$, so the service times add little noise), but they come for free.

7.5. Stratified Sampling

Goal: estimate $\theta = \mathbf{E}[X]$, using an intermediate r.v. $Y \in \{1, 2, \dots, k\}$ with known p.m.f.: $P(Y = i)$, $i = 1, \dots, k$, is given. By the law of total expectation,

$$\theta = \mathbf{E}[X] = \sum_{i=1}^k \mathbf{E}[X \mid Y = i]\, P(Y = i) .$$

Let $\bar X_i$ be an estimator of $\mathbf{E}[X \mid Y = i]$: the average of $n_i$ observations of $X$ generated conditional on $Y = i$ (from “stratum” $i$). Then the stratified sampling estimator is

$$\hat\theta = \sum_{i=1}^k \bar X_i\, P(Y = i) .$$

Why it helps: the strata are sampled independently, so

$$\operatorname{Var}(\hat\theta) = \sum_{i=1}^k P(Y = i)^2\,\frac{\operatorname{Var}(X \mid Y = i)}{n_i} .$$

With $N$ runs in total and proportional allocation, $n_i = N\,P(Y = i)$, this is $\frac{1}{N}\sum_i P(Y = i)\operatorname{Var}(X \mid Y = i) = \frac{1}{N}\mathbf{E}[\operatorname{Var}(X \mid Y)] \le \frac{1}{N}\operatorname{Var}(X)$, the crude Monte Carlo variance with $N$ runs. Stratification removes the part $\frac{1}{N}\operatorname{Var}(\mathbf{E}[X \mid Y])$ of the variance, i.e. the randomness of how many observations happen to fall in each stratum (law of total variance, §7.6).

Ex: Estimate $\mathbf{E}[g(U)] = \int_0^1 g(x)\, dx$, $U \sim U(0, 1)$. Two approaches, each with $2n$ evaluations of $g$.

Crude Monte Carlo based on a sample of $2n$ runs:

$$\hat g = \frac{1}{2n}\sum_{i=1}^{2n} g(U_i),$$

$$\operatorname{Var}(\hat g) = \frac{1}{4n^2}\sum_{i=1}^{2n}\operatorname{Var}(g(U_i)) = \frac{1}{2n}\operatorname{Var}(g(U)) = \frac{1}{2n}\left[\int_0^1 g^2(x)\, dx - \left(\int_0^1 g(x)\, dx\right)^2\right].$$

Stratified sampling with

$$Y = \begin{cases} 0, & U \in (0, \frac{1}{2}), \\ 1, & U \in (\frac{1}{2}, 1), \end{cases} \qquad P(Y = 0) = P(Y = 1) = \tfrac{1}{2} .$$

$$\begin{aligned} \mathbf{E}[g(U)] &= \mathbf{E}\big[g(U) \mid U \in (0, \tfrac{1}{2})\big]\cdot\tfrac{1}{2} + \mathbf{E}\big[g(U) \mid U \in (\tfrac{1}{2}, 1)\big]\cdot\tfrac{1}{2} \\ &= \mathbf{E}[g(\bar U)]\cdot\tfrac{1}{2} + \mathbf{E}[g(\tilde U)]\cdot\tfrac{1}{2}, \qquad \bar U \sim U(0, \tfrac{1}{2}),\ \tilde U \sim U(\tfrac{1}{2}, 1) \\ &= \mathbf{E}\big[g(\tfrac{1}{2}U)\big]\cdot\tfrac{1}{2} + \mathbf{E}\big[g(\tfrac{1}{2} + \tfrac{1}{2}U)\big]\cdot\tfrac{1}{2}, \end{aligned}$$

since $\frac{1}{2}U \sim U(0, \frac{1}{2})$ and $\frac{1}{2} + \frac{1}{2}U \sim U(\frac{1}{2}, 1)$. Estimator ($n$ points in each half):

$$\tilde g = \frac{1}{2}\cdot\frac{1}{n}\sum_{i=1}^{n} g\big(\tfrac{1}{2}U_i\big) + \frac{1}{2}\cdot\frac{1}{n}\sum_{i=n+1}^{2n} g\big(\tfrac{1}{2} + \tfrac{1}{2}U_i\big),$$

$$\operatorname{Var}(\tilde g) = \operatorname{Var}(\hat g) - \frac{1}{2n}\left[\int_0^{1/2} g(x)\, dx - \int_{1/2}^1 g(x)\, dx\right]^2 \le \operatorname{Var}(\hat g) .$$

Derivation: let $A = \int_0^{1/2} g$, $B = \int_{1/2}^1 g$, $A_2 = \int_0^{1/2} g^2$, $B_2 = \int_{1/2}^1 g^2$. The density of $\frac{1}{2}U$ is 2 on $(0, \frac{1}{2})$, so $\mathbf{E}[g(\frac{1}{2}U)] = 2A$, $\mathbf{E}[g^2(\frac{1}{2}U)] = 2A_2$, and $\operatorname{Var}(g(\frac{1}{2}U)) = 2A_2 - 4A^2$; likewise for the other half. Then

$$\operatorname{Var}(\tilde g) = \frac{1}{4n}\big(2A_2 - 4A^2\big) + \frac{1}{4n}\big(2B_2 - 4B^2\big) = \frac{1}{2n}\big[A_2 + B_2 - 2A^2 - 2B^2\big],$$

$$\operatorname{Var}(\hat g) = \frac{1}{2n}\big[A_2 + B_2 - (A + B)^2\big],$$

$$\operatorname{Var}(\hat g) - \operatorname{Var}(\tilde g) = \frac{1}{2n}\big[2A^2 + 2B^2 - (A + B)^2\big] = \frac{1}{2n}(A - B)^2 \ge 0 .$$

The reduction is large when the two halves have very different integrals. Crude Monte Carlo may by chance put more points in one half; stratification puts exactly $n$ in each:

E.g. $g(x) = e^x$: $A = e^{1/2} - 1 = 0.649$, $B = e - e^{1/2} = 1.070$, $(A - B)^2 = 0.177$, and $\operatorname{Var}(g(U)) = 0.242$. So $\operatorname{Var}(\hat g) = 0.242/2n$ and $\operatorname{Var}(\tilde g) = (0.242 - 0.177)/2n = 0.065/2n$: a 73% reduction.

7.6. Conditioning

Goal: estimate $\mu = \mathbf{E}[X]$, using an intermediate r.v. $Z$ satisfying: if the value of $Z$ is known, then $\mathbf{E}[X \mid Z = z]$ is a known, deterministic function of $z$ (it can be computed exactly, without simulation).

$$\mu = \mathbf{E}[X] = \mathbf{E}\big[\mathbf{E}[X \mid Z]\big] .$$

The conditioning estimator is

$$\hat\mu = \mathbf{E}[X \mid Z = z],$$

evaluated at the observed value $z$ of $Z$ (i.e. the random variable $\mathbf{E}[X \mid Z]$): instead of using the observed $X$, observe $Z$ and use the exact conditional mean.

Ex: a single-server queue (FCFS) with exponential service times. $\mu = \mathbf{E}[D_i]$, $D_i$ = delay of the $i$-th customer. $Z$ = number of customers in the system seen by the $i$-th arrival. Then

$$\mathbf{E}[D_i \mid Z = z] = z\,\mathbf{E}[S] .$$

Why: customer $i$ has to wait for the $z$ customers ahead of him. The $z - 1$ customers waiting in the queue each need a full service time (mean $\mathbf{E}[S]$). The customer in service needs its remaining service time, which by the memoryless property of the exponential distribution also has mean $\mathbf{E}[S]$, no matter how long the service has already lasted. In total, $z\,\mathbf{E}[S]$. (With non-exponential service times the remaining service time would depend on the elapsed service time, and the formula would fail.) So the conditioning estimator of $\mathbf{E}[D_i]$ is $Z\,\mathbf{E}[S]$, and for the average delay of the first $n$ customers it is $\frac{1}{n}\sum_{i=1}^n Z_i\,\mathbf{E}[S]$.

Law of total variance:

$$\operatorname{Var}(X) = \mathbf{E}\big[\operatorname{Var}(X \mid Z)\big] + \operatorname{Var}\big(\mathbf{E}[X \mid Z]\big)$$

(total variance = average variance within each value of $Z$ + variance of the conditional means). Hence

$$\operatorname{Var}\big(\mathbf{E}[X \mid Z]\big) = \operatorname{Var}(X) - \mathbf{E}\big[\operatorname{Var}(X \mid Z)\big] \le \operatorname{Var}(X),$$

since $\mathbf{E}[\operatorname{Var}(X \mid Z)] \ge 0$. So the conditioning estimator $\hat\mu = \mathbf{E}[X \mid Z]$ is unbiased (because $\mu = \mathbf{E}[X] = \mathbf{E}[\mathbf{E}[X \mid Z]]$) and its variance is never larger than that of $X$. The randomness it removes is $\mathbf{E}[\operatorname{Var}(X \mid Z)]$: here, the randomness of the actual service times of the customers ahead, which is averaged out exactly instead of being simulated.

E.g. M/M/1 queue ($\lambda = 1$, $\mathbf{E}[S] = 0.9$, empty start), $X$ = average delay of the first 100 customers (20,000 runs): both estimators have mean 4.13, but $\operatorname{Var}(X) = 8.48$ while $\operatorname{Var}\big(\frac{1}{100}\sum_i 0.9\,Z_i\big) = 6.56$, a 23% reduction.

7.7. Importance Sampling

Goal: estimate

$$\mu = \mathbf{E}[X] = \mathbf{E}[h(Y)] = \int_{-\infty}^{\infty} h(y) f(y)\, dy = \int_{-\infty}^{\infty} h(y)\,\frac{f(y)}{g(y)}\, g(y)\, dy,$$

where $Y$ has density $f$, and $g$ is another density that we choose (and can sample from). We need $f$ to be absolutely continuous with respect to $g$ (the board said “$g$ with respect to $f$”, but the “i.e.” part is the right condition), i.e. whenever $g(y) = 0$, $f(y) = 0$. Equivalently:

$$\text{support of } f \subseteq \text{support of } g,$$

where the support of $f$ is the set of values of $y$ s.t. $f(y) > 0$.
Why this is needed: we divide by $g(y)$, and any region where $f(y) > 0$ but $g(y) = 0$ would never be sampled from $g$, so that part of the integral would be missed. (Strictly, it is enough where $h(y) \ne 0$.)

Write $\mathbf{E}_f$ for an expectation when $Y$ has density $f$, and $\mathbf{E}_g$ when $Y$ has density $g$. Then

$$\mu = \mathbf{E}[X] = \mathbf{E}_g\!\left[h(Y)\,\frac{f(Y)}{g(Y)}\right].$$

Generate $Y_1, Y_2, \dots, Y_n \sim g$ (i.i.d.) and use

$$\hat\mu(n) = \frac{1}{n}\sum_{i=1}^n h(Y_i)\underbrace{\left\{\frac{f(Y_i)}{g(Y_i)}\right\}}_{\text{likelihood ratio}} .$$

This one is unbiased, $\mathbf{E}_g[\hat\mu(n)] = \mu$, just like the crude estimator $\frac{1}{n}\sum_{i=1}^n h(Y_i)$ with $Y_i \sim f$. The likelihood ratio is a weight that corrects for sampling from the “wrong” density: values that $g$ produces more often than $f$ would ($g > f$) get weight $< 1$, values that $g$ produces less often get weight $> 1$.

Variance:

$$\begin{aligned} \operatorname{Var}(\hat\mu(n)) &= \frac{1}{n}\left(\mathbf{E}_g\!\left[\frac{h^2(Y)f^2(Y)}{g^2(Y)}\right] - \underbrace{\left\{\mathbf{E}_g\!\left[\frac{h(Y)f(Y)}{g(Y)}\right]\right\}^2}_{=\ \mu^2}\right) \\ &= \frac{1}{n}\,\mathbf{E}_g\!\left[\frac{h^2(Y)f^2(Y)}{g^2(Y)}\right] - \frac{1}{n}\mu^2 \\ &= \frac{1}{n}\int \frac{h^2(y)f^2(y)}{g^2(y)}\, g(y)\, dy - \frac{1}{n}\mu^2 \\ &= \frac{1}{n}\int \frac{h^2(y)f^2(y)}{g(y)}\, dy - \frac{1}{n}\mu^2 \qquad (\text{one } g(y) \text{ cancels}) \\ &= \frac{1}{n}\int h^2(y)\,\frac{f(y)}{g(y)}\, f(y)\, dy - \frac{1}{n}\mu^2 \\ &= \frac{1}{n}\,\mathbf{E}_f\!\left[\frac{h^2(Y)f(Y)}{g(Y)}\right] - \frac{1}{n}\mu^2 . \end{aligned}$$

“$Y \sim f$? Wait, but $Y$ is already $\sim g$?” In the simulation, the $Y_i$’s are always generated from $g$; nothing changes there. The subscript of $\mathbf{E}$ only says which density the integral is taken against: $\mathbf{E}_g[k(Y)] = \int k(y) g(y)\, dy$ and $\mathbf{E}_f[k(Y)] = \int k(y) f(y)\, dy$, and the $Y$ inside is just a name for the integration variable. The last step does not change the number: it regroups the same integral $\int h^2 f^2/g\, dy$, using $f(y)$ instead of $g(y)$ as the density, so that it can be compared with crude Monte Carlo:

$$\text{crude } (Y_i \sim f):\ \ \operatorname{Var} = \frac{1}{n}\,\mathbf{E}_f\big[h^2(Y)\big] - \frac{1}{n}\mu^2,$$

$$\text{importance sampling } (Y_i \sim g):\ \ \operatorname{Var} = \frac{1}{n}\,\mathbf{E}_f\!\left[h^2(Y)\,\frac{f(Y)}{g(Y)}\right] - \frac{1}{n}\mu^2 .$$

The $\mu^2$ terms are the same, so importance sampling wins when $f(y)/g(y) < 1$ mostly where $h^2(y) f(y)$ is large: $g$ should generate the “important” values (those that contribute most to the integral) more often than $f$ does. Hence the name.

Choosing $g$: only the first term depends on $g$, so we want

$$\min_g\ \mathbf{E}_f\!\left[\frac{h^2(Y)f(Y)}{g(Y)}\right].$$

• The best possible choice (for $h \ge 0$) is $g^*(y) = h(y)f(y)/\mu$: then $h(y)f(y)/g^*(y) = \mu$ for every $y$ and the variance is 0. It cannot be used, because it requires the unknown $\mu$, but it is the guideline: choose $g$ roughly proportional to $h \cdot f$.
• In practice, $g$ is chosen from a family and its parameter is optimized. E.g. if $f$ is $\text{Exp}(\lambda)$, take $g = \text{Exp}(\lambda')$ and choose the rate $\lambda'$. (The board called this rate $\mu$; I renamed it to avoid a clash with $\mu = \mathbf{E}[X]$.)

Ex (worked out for this family): estimate $p = P(Y > a)$ for $Y \sim \text{Exp}(\lambda)$, so $h(y) = I\{y > a\}$ (the indicator $I$ is explained in the next example). For $\lambda' < 2\lambda$,

$$\mathbf{E}_f\!\left[\frac{h^2(Y)f(Y)}{g(Y)}\right] = \int_a^\infty \frac{(\lambda e^{-\lambda y})^2}{\lambda' e^{-\lambda' y}}\, dy = \frac{\lambda^2\, e^{-(2\lambda - \lambda')a}}{\lambda'(2\lambda - \lambda')},$$

and setting the derivative of its logarithm to 0 gives

$$\lambda'^* = \lambda + \frac{1}{a} - \sqrt{\lambda^2 + \frac{1}{a^2}} \;\approx\; \frac{1}{a} \quad \text{for large } a,$$

i.e. the best $g$ has mean $\approx a$: it puts the samples around the threshold. With $\lambda = 1$ and $a = 20$, $p = e^{-20} = 2.06 \times 10^{-9}$. One crude sample has standard deviation $\sqrt{p(1-p)} \approx 22{,}000 \times p$. One importance sample with $\lambda'^* = 0.049$ has standard deviation $5.2 \times p$.

Ex: $X \sim N(0, 1)$. Estimate $p = P(10 < X < 11)$ by simulation. The true value is

$$p = P(10 < X < 11) = 7.6 \times 10^{-24} \quad (\text{board: } \approx 10^{-24}, \text{ the right order of magnitude}).$$

Monte Carlo simulation: generate $X_1, \dots, X_n \sim N(0, 1)$. Write the probability as an expectation:

$$P(10 < X < 11) = \mathbf{E}\big[I\{10 < X < 11\}\big] .$$

What is $I$? The indicator function of an event: $I\{A\} = 1$ if $A$ happens and $0$ otherwise. So $I\{10 < X < 11\}$ is a Bernoulli r.v.: it equals 1 with probability $p$ and 0 with probability $1 - p$. Hence $\mathbf{E}[I] = 1 \cdot p + 0 \cdot (1 - p) = p$ and $\operatorname{Var}(I) = p(1 - p)$. This turns any probability into an expectation, with $h(y) = I\{10 < y < 11\}$. The estimator

$$\tilde\mu(n) = \frac{1}{n}\sum_{i=1}^n I\{10 < X_i < 11\}$$

is the fraction of the $X_i$’s that fall in $(10, 11)$, and

$$\operatorname{Var}(\tilde\mu(n)) = \frac{p(1 - p)}{n} \approx \frac{10^{-24}}{n} .$$

A 95% CI for $\mu = p$ is

$$\left[\tilde\mu(n) - 1.96\sqrt{\frac{10^{-24}}{n}},\;\; \tilde\mu(n) + 1.96\sqrt{\frac{10^{-24}}{n}}\right].$$

• $n = 10^6$: half-width $= 1.96\sqrt{10^{-24}/10^{6}} = 1.96 \times 10^{-15}$, about $10^9$ times larger than $p$.
• $n = 10^{10}$: half-width $= 1.96\sqrt{10^{-24}/10^{10}} = 1.96\sqrt{10^{-34}} = 1.96 \times 10^{-17}$. This is the $10^{-17}$ on the board: the half-width is still about ten million times larger than the true mean.

“Variance is too large compared with the true probability?” Yes, exactly. The variance $10^{-34}$ is tiny in absolute terms, but what matters is the error relative to the tiny $p$:

$$\frac{\text{standard deviation of } \tilde\mu(n)}{p} = \frac{\sqrt{p(1-p)/n}}{p} \approx \frac{1}{\sqrt{np}} .$$

For the half-width to be even as small as $p$ itself (a 100% error!), we need $1.96/\sqrt{np} \le 1$, i.e. $n \ge 1.96^2/p \approx 5 \times 10^{23}$ runs. Intuitively: with $n = 10^{10}$ runs, the expected number of $X_i$’s that land in $(10, 11)$ is $np \approx 10^{-13}$. In practice no $X_i$ ever lands there, so $\tilde\mu(n) = 0$, the sample variance is 0 too, and the computed CI is $[0, 0]$: a meaningless answer. Crude Monte Carlo cannot see an event that it never generates.

Importance sampling. Sample from $g = N(10, 1)$ instead, i.e. move the normal density to the interval:

Write $\phi$ for the $N(0, 1)$ density ($= f$) and $\phi_{10}$ for the $N(10, 1)$ density ($= g$); “$N(0,1)/N(10,1)$” on the board means the ratio $\phi(x)/\phi_{10}(x)$ of the two densities. Then

$$\begin{aligned} \mu = P(10 < X < 11) &= \mathbf{E}\big[I\{10 < X < 11\}\big] \\ &= \int_{-\infty}^{\infty} I\{10 < x < 11\}\,\frac{\phi(x)}{\phi_{10}(x)}\,\underbrace{\phi_{10}(x)}_{g(x)}\, dx \\ &= \mathbf{E}_g\!\left[I\{10 < X < 11\}\,\frac{\phi(X)}{\phi_{10}(X)}\right] \\ &= \mathbf{E}_g\!\left[I\{10 < X < 11\}\,\frac{\frac{1}{\sqrt{2\pi}}\,e^{-X^2/2}}{\frac{1}{\sqrt{2\pi}}\,e^{-(X-10)^2/2}}\right] \\ &= \mathbf{E}_g\!\left[I\{10 < X < 11\}\, e^{-X^2/2 + (X-10)^2/2}\right] \\ &= \mathbf{E}_g\!\left[I\{10 < X < 11\}\, e^{-10X + 50}\right]. \end{aligned}$$

The $\frac{1}{\sqrt{2\pi}}$'s cancel out, and $-\frac{x^2}{2} + \frac{(x - 10)^2}{2} = -\frac{x^2}{2} + \frac{x^2 - 20x + 100}{2} = -10x + 50$.

Generate $X_1, \dots, X_n \sim N(10, 1)$. (The board says $N(0, 1)$, but the samples must come from $g = N(10, 1)$: that is the whole point.) The estimator is

$$\hat\mu(n) = \frac{1}{n}\sum_{i=1}^n I\{10 < X_i < 11\}\, e^{-10X_i + 50}$$

(written $\tilde\mu(n)$ on the board; I use $\hat\mu(n)$, as in the general formula above, to keep it apart from the crude estimator $\tilde\mu(n)$).

Variance:

$$\begin{aligned} \operatorname{Var}(\hat\mu(n)) &= \frac{1}{n}\Big(\mathbf{E}_g\big[I^2\{10 < X < 11\}\, e^{-20X + 100}\big] - p^2\Big) \\ &= \frac{1}{n}\Big(\mathbf{E}_g\big[I\{10 < X < 11\}\, e^{-20X + 100}\big] - p^2\Big) \\ &\le \frac{1}{n}\, e^{-100} - \frac{1}{n}\, p^2 = \frac{1}{n}\big(10^{-44} - 10^{-48}\big) \approx \frac{1}{n}\, 10^{-44} . \end{aligned}$$

Step by step:

1. This is the general formula $\frac{1}{n}\big(\mathbf{E}_g[(\text{weighted observation})^2] - \mu^2\big)$, and the square of the weight is $(e^{-10x + 50})^2 = e^{-20x + 100}$.
2. $I^2 = I$, because $0^2 = 0$ and $1^2 = 1$.
3. The bound: the indicator is nonzero only when $10 < X < 11$, and there $-20X + 100 < -200 + 100 = -100$. So $I\{10 < X < 11\}\, e^{-20X + 100} \le e^{-100}$ always, and its expectation is $\le e^{-100} = 3.7 \times 10^{-44} \approx 10^{-44}$.
4. $p^2 \approx (10^{-24})^2 = 10^{-48}$ (with the exact $p$, $5.8 \times 10^{-47}$), negligible.

Compare with crude Monte Carlo:

$$\operatorname{Var}(\tilde\mu(n)) \approx \frac{10^{-24}}{n} \qquad\text{vs.}\qquad \operatorname{Var}(\hat\mu(n)) \lesssim \frac{10^{-44}}{n} :$$

the variance is about $10^{20}$ times smaller, for the same $n$. In relative terms, the standard deviation of one sample is now at most $\sqrt{10^{-44}}/10^{-24} = 100$ times $p$ (bound with $p \approx 10^{-24}$), instead of $1/\sqrt p \approx 10^{12}$ times $p$. E.g. with $n = 10^6$ the 95% half-width is at most $1.96\sqrt{10^{-44}/10^{6}} \approx 2 \times 10^{-25}$, a fraction of $p$, while crude Monte Carlo gave $2 \times 10^{-15}$ with the same $n$.

The bound is conservative. The exact per-sample variance is $6.8 \times 10^{-46}$ (computed numerically), i.e. a standard deviation of $3.4 \times p$. E.g. with $n = 10^4$: 3421 of the $X_i$’s land in $(10, 11)$, $\hat\mu = 7.44 \times 10^{-24}$, half-width $5.0 \times 10^{-25}$, so the 95% CI is $[6.94,\ 7.95] \times 10^{-24}$, which covers $p = 7.62 \times 10^{-24}$. Over 2000 repetitions, this CI covered $p$ in 94.3% of cases.

Why it works: under $g$, about 34% of the $X_i$’s land in $(10, 11)$ ($P(10 < X < 11) = \Phi(1) - \Phi(0) = 0.341$ for $X \sim N(10, 1)$), instead of $np \approx 10^{-13}$ expected hits with crude sampling. Each hit is then multiplied by the likelihood ratio $e^{-10x + 50}$, between $e^{-60} = 8.8 \times 10^{-27}$ (at $x = 11$) and $e^{-50} = 1.9 \times 10^{-22}$ (at $x = 10$), which scales the frequent hits back down to the true, tiny probability.

Choice of the shift: among $g = N(\theta, 1)$, the variance is smallest at $\theta \approx 10.05$ ($3.43 \times p$), so $\theta = 10$ is essentially optimal ($\theta = 10.5$ gives $3.6 \times p$, $\theta = 11$ gives $4.4 \times p$). The reason: $f$ falls so fast that almost all of the probability of $(10, 11)$ sits just above 10, i.e. $h \cdot f$ is concentrated near 10, and the guideline “$g$ roughly proportional to $h \cdot f$” says to put the samples there.