25. Mean of a Likelihood Ratio Process#
GPU
This lecture was built using a machine with access to a GPU — although it will also run without one.
Google Colab has a free tier with GPUs that you can access as follows:
Click on the “play” icon top right
Select Colab
Set the runtime environment to include a GPU
25.1. Overview#
In Likelihood Ratio Processes we described a peculiar property of a likelihood ratio process, namely, that its mean equals one for all \(t \geq 0\) despite its converging to zero almost surely.
While it is easy to verify that peculiar property analytically (i.e., in population), it is challenging to use a computer simulation to verify it via an application of a law of large numbers that entails studying sample averages of repeated simulations.
To confront this challenge, this lecture puts importance sampling to work to accelerate convergence of sample averages to population means.
We use importance sampling to estimate the mean of a cumulative likelihood ratio \(L\left(\omega^t\right) = \prod_{i=1}^t \ell \left(\omega_i\right)\).
In addition to what’s in Anaconda, this lecture will need the following libraries:
!pip install jax
We start by importing some Python packages.
import jax
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt
from jax.scipy.stats import beta
from typing import NamedTuple
from functools import partial
# Set JAX to use 64-bit floats
jax.config.update("jax_enable_x64", True)
25.2. Mathematical expectation of likelihood ratio#
In Likelihood Ratio Processes, we studied a likelihood ratio \(\ell \left(\omega_t\right)\)
where \(f\) is the Beta\((F_a, F_b) = \) Beta\((1, 1)\) density, the uniform density on \([0,1]\), and \(g\) is the Beta\((G_a, G_b) = \) Beta\((3, 1.2)\) density.
Assume that an i.i.d. random variable \(\omega_t \in \Omega\) is generated by \(g\).
The cumulative likelihood ratio \(L \left(\omega^t\right)\) is
Our goal is to approximate the mathematical expectation \(E \left[ L\left(\omega^t\right) \right]\) well.
In this section of Likelihood Ratio Processes, we showed that \(E \left[ L\left(\omega^t\right) \right]\) equals \(1\) for all \(t\).
We want to check out how well this holds if we replace \(E\) with sample averages from simulations.
This turns out to be easier said than done because for Beta distributions assumed above, \(L\left(\omega^t\right)\) has a very skewed distribution with a very long tail as \(t \rightarrow \infty\).
This property makes it difficult efficiently and accurately to estimate the mean by standard Monte Carlo simulation methods.
In this lecture we explore how a standard Monte Carlo method fails.
We also show how importance sampling provides a more computationally efficient way to approximate the mean of the cumulative likelihood ratio.
We first take a look at the density functions f and g .
# Parameters for the model
class ImpSampleParams(NamedTuple):
F_a: float = 1.0 # Beta parameters for f
F_b: float = 1.0
G_a: float = 3.0 # Beta parameters for g
G_b: float = 1.2
params = ImpSampleParams()
def f(w):
return beta.pdf(w, params.F_a, params.F_b)
def g(w):
return beta.pdf(w, params.G_a, params.G_b)
w_range = np.linspace(1e-2, 1-1e-5, 1000)
plt.plot(w_range, g(w_range), lw=2, label='g')
plt.plot(w_range, f(w_range), lw=2, label='f')
plt.xlabel(r'$\omega$')
plt.legend()
plt.show()
Fig. 25.1 Beta density functions \(f\) and \(g\)#
The likelihood ratio is l(w)=f(w)/g(w).
def l(w):
return f(w) / g(w)
plt.plot(w_range, l(w_range), lw=2)
plt.yscale('log')
plt.xlabel(r'$\omega$')
plt.ylabel(r'$\ell(\omega)$')
plt.show()
Fig. 25.2 Likelihood ratio \(\ell(\omega)\), log scale#
Note the log scale on the vertical axis: on a linear scale the entire curve would be flattened against zero by the spike at the left-hand edge.
Fig. 25.1 shows that as \(\omega \rightarrow 0\), \(f \left(\omega\right)\) is unchanged while \(g \left(\omega\right) \rightarrow 0\).
Hence the likelihood ratio in Fig. 25.2 diverges to infinity, at rate \(\omega^{-2}\).
The same thing happens as \(\omega \rightarrow 1\), since \(G_b > 1\) also forces \(g\left(\omega\right) \rightarrow 0\) there, which explains the upturn at the right-hand edge.
But that divergence is far slower — at rate \(\left(1-\omega\right)^{-1/5}\) — and it is the behavior near \(\omega = 0\) that causes the trouble below.
A Monte Carlo approximation of \(E \left[L\left(\omega^t\right)\right]\) would repeatedly draw a sequence \(\omega^t = \left(\omega_1, \ldots, \omega_t\right)\) of \(t\) independent observations from \(g\), form the product \(L\left(\omega^t\right) = \prod_{i=1}^t \ell \left(\omega_i\right)\) for each such sequence, then average these products across independently drawn sequences.
Because \(g(\omega) \rightarrow 0\) as \(\omega \rightarrow 0\), such a simulation procedure undersamples a part of the sample space \([0,1]\) that it is important to visit often in order to do a good job of approximating the mathematical expectation of \(\ell \left(\omega\right)\).
Every one of the \(t\) factors in \(L\left(\omega^t\right)\) is distorted in the same way, so the problem compounds as \(t\) grows.
We illustrate this numerically below.
25.3. Importance sampling#
We circumvent the issue by using a change of distribution called importance sampling.
Instead of drawing from \(g\) to generate data during the simulation, we use an alternative distribution \(h\) to generate draws of \(\omega\).
In Likelihood Ratio Processes, \(h\) denotes a third density that generates the data; here \(h\) denotes an importance distribution that we choose.
The idea is to design \(h\) so that it oversamples the region of \(\Omega\) where \(\ell \left(\omega_t\right)\) has large values but low density under \(g\).
After we construct a sample in this way, we must then weight each realization by the likelihood ratio of \(g\) and \(h\) when we compute the empirical mean of the likelihood ratio.
By doing this, we properly account for the fact that we are using \(h\) and not \(g\) to simulate data.
To illustrate, suppose we were interested in \({E}\left[\ell\left(\omega\right)\right]\).
We could simply compute:
where \(\omega_i^g\) indicates that \(\omega_i\) is drawn from \(g\).
But using our insight from importance sampling, we could instead calculate the object:
where \(\omega_i\) is now drawn from importance distribution \(h\).
Notice that the above two are exactly the same population objects:
The second equality requires \(h(\omega) > 0\) wherever \(\ell(\omega) g(\omega) = f(\omega) > 0\).
Since we must use an \(h\) that has larger mass in parts of the distribution to which \(g\) puts low mass, we use \(h=Beta(0.5, 0.5)\) as our importance distribution.
The plots compare \(g\) and \(h\).
g_a, g_b = params.G_a, params.G_b
h_a, h_b = 0.5, 0.5
w_range = np.linspace(1e-5, 1-1e-5, 1000)
plt.plot(w_range, g(w_range),
lw=2, label=f'g=Beta({g_a}, {g_b})')
plt.plot(w_range, beta.pdf(w_range, 0.5, 0.5),
lw=2, label=f'h=Beta({h_a}, {h_b})')
plt.legend()
plt.ylim([0., 3.])
plt.show()
Fig. 25.3 Data and importance sampling distributions#
25.4. Approximating a cumulative likelihood ratio#
We now study how to use importance sampling to approximate \({E} \left[L\left(\omega^t\right)\right] = E \left[\prod_{i=1}^t \ell \left(\omega_i\right)\right]\).
To write code for general Beta densities, we let \(p\) denote the density under which we want the expectation and \(q\) the density from which we draw.
In our application \(p = g\) and \(q\) is an importance distribution such as \(h\).
In Likelihood Ratio Processes, \(q\) denotes the density that nature chose; here \(q\) denotes the importance distribution.
As above, our plan is to draw sequences \(\omega^t\) from \(q\) and then re-weight the likelihood ratio appropriately.
The change of distribution holds in population, exactly as it did for a single draw:
This suggests the estimator
where \(\omega_{n,i}^q\) is the \(i\)-th observation of the \(n\)-th sequence, drawn from the importance distribution \(q\).
Here \(\frac{p\left(\omega_{n,i}^q\right)}{q\left(\omega_{n,i}^q\right)}\) is the weight we assign to each data point \(\omega_{n,i}^q\).
Below we prepare a Python function for computing the importance sampling estimates given any beta distributions \(p\), \(q\).
def estimate_single_path(key, p_a, p_b, q_a, q_b, T):
"""
Estimation for a single sample path.
"""
def loop_body(i, carry):
L, weight, key_state = carry
key_state, subkey = jax.random.split(key_state)
w = jax.random.beta(subkey, q_a, q_b)
# Keep draws off the boundary, where the log densities
# below produce the undefined form 0 * inf
w = jnp.clip(w, 1e-12, 1 - 1e-12)
# Compute likelihood ratio using f/g functions
likelihood_ratio = f(w) / g(w)
L = L * likelihood_ratio
# Importance sampling weight
p_w = beta.pdf(w, p_a, p_b)
q_w = beta.pdf(w, q_a, q_b)
weight = weight * (p_w / q_w)
return (L, weight, key_state)
# Use fori_loop for dynamic T values
final_L, final_weight, _ = jax.lax.fori_loop(
0, T, loop_body, (1.0, 1.0, key)
)
return final_L * final_weight
@partial(jax.jit, static_argnames=['N'])
def estimate(key, p_a, p_b, q_a, q_b, T=1, N=10000):
"""Estimation of a batch of sample paths."""
keys = jax.random.split(key, N)
# Vectorize over keys, holding the parameters fixed
estimates = jax.vmap(
estimate_single_path,
in_axes=(0, None, None, None, None, None)
)(keys, p_a, p_b, q_a, q_b, T)
return jnp.mean(estimates)
Consider the case when \(T=1\), which amounts to approximating \(E\left[\ell\left(\omega\right)\right]\).
For the standard Monte Carlo estimate, we can set \(p=g\) and \(q=g\).
estimate(jax.random.key(0), g_a, g_b, g_a, g_b,
T=1, N=10000)
Array(0.98300826, dtype=float64)
For our importance sampling estimate, we set \(q = h\).
estimate(jax.random.key(1), g_a, g_b, h_a, h_b,
T=1, N=10000)
Array(1.00722995, dtype=float64)
Evidently, even at \(T=1\), our importance sampling estimate is closer to \(1\) than is the Monte Carlo estimate.
Bigger differences arise when computing expectations over longer sequences, \(E\left[L\left(\omega^t\right)\right]\).
Setting \(T=10\), we find that the Monte Carlo method severely underestimates the mean while importance sampling still produces an estimate close to its theoretical value of unity.
estimate(jax.random.key(2), g_a, g_b, g_a, g_b,
T=10, N=10000)
Array(0.60388316, dtype=float64)
estimate(jax.random.key(3), g_a, g_b, h_a, h_b,
T=10, N=10000)
Array(0.99185478, dtype=float64)
The Monte Carlo method underestimates because the likelihood ratio \(L(\omega^T) = \prod_{t=1}^T \frac{f(\omega_t)}{g(\omega_t)}\) has a highly skewed distribution under \(g\).
Most samples from \(g\) produce small likelihood ratios, while the true mean requires occasional very large values that are rarely sampled.
In our case, since \(g(\omega) \to 0\) as \(\omega \to 0\) while \(f(\omega)\) remains constant, the Monte Carlo procedure undersamples precisely where the likelihood ratio \(\frac{f(\omega)}{g(\omega)}\) is largest.
In fact the situation is worse than skewness — the Monte Carlo estimator has infinite variance.
To see this, note that the second moment of a single likelihood ratio is
Here \(f\) is the uniform density, while \(g\left(\omega\right)\) vanishes like \(\omega^{G_a - 1} = \omega^2\) as \(\omega \to 0\).
The integrand therefore behaves like \(\omega^{-2}\) near the origin, and the integral diverges.
So \(\ell\left(\omega\right)\) has no finite variance, and since the draws are independent, neither does \(L\left(\omega^t\right)\) for any \(t\).
Likelihood Ratio Processes meets the same divergent integral as the mean of \(\ell(\omega)\) under \(f\), since \(E^g\left[\ell(\omega)^2\right] = E^f\left[\ell(\omega)\right]\).
The mean \(E^g\left[\ell(\omega)\right] = 1\) is finite, so the strong law of large numbers still applies: the Monte Carlo average converges to \(1\) almost surely as the number of draws \(N\) grows.
The ordinary central limit theorem does not apply, and its \(N^{-1/2}\) rate of convergence fails.
Near the origin, \(g(\omega) \approx k \omega^2\) with \(k = 1/B(3, 1.2)\), where \(B\) is the Beta function.
For large \(x\), the event \(\ell(\omega) > x\) occurs mainly when \(\omega < (k x)^{-1/2}\), an event with probability about \(\frac{1}{3} k^{-1/2} x^{-3/2}\) under \(g\).
A tail that decays like \(x^{-3/2}\) puts \(\ell(\omega)\) in the domain of attraction of a stable law with index \(3/2\).
By the generalized central limit theorem, the error of a sample mean of \(N\) draws shrinks like \(N^{-1/3}\), and its limiting distribution is skewed, with a long right tail.
Exercise 25.2 confirms the tail index and the \(N^{-1/3}\) rate numerically.
As \(T\) increases, the problem worsens, making standard Monte Carlo increasingly unreliable.
Importance sampling with \(q = h\) fixes the infinite variance by sampling more often near \(\omega = 0\), where \(\ell\) is large.
25.5. Distribution of sample mean#
We next study the bias and efficiency of the Monte Carlo and importance sampling approaches.
The code below repeats the estimate N_simu times, so that we can look at the distribution of the estimates that each method produces.
@partial(jax.jit, static_argnames=['N_simu', 'N_samples'])
def simulate(key, p_a, p_b, q_a, q_b, N_simu, T=1,
N_samples=10000):
"""Repeat the estimate N_simu times, drawing from q."""
keys = jax.random.split(key, N_simu)
return jax.vmap(
lambda k: estimate(k, p_a, p_b, q_a, q_b, T,
N_samples)
)(keys)
Setting \(q = p\) recovers standard Monte Carlo, since every importance weight is then identically one.
Again, we first consider estimating \({E} \left[\ell\left(\omega\right)\right]\) by setting T=1.
We simulate \(1000\) times for each method.
N_simu = 1000
μ_L_g = simulate(jax.random.key(4), g_a, g_b,
g_a, g_b, N_simu)
μ_L_h = simulate(jax.random.key(5), g_a, g_b,
h_a, h_b, N_simu)
# standard Monte Carlo (mean and variance)
jnp.mean(μ_L_g), jnp.var(μ_L_g)
(Array(0.99558545, dtype=float64), Array(0.00448641, dtype=float64))
# importance sampling (mean and variance)
jnp.mean(μ_L_h), jnp.var(μ_L_h)
(Array(0.99985196, dtype=float64), Array(2.46227759e-05, dtype=float64))
Although both methods tend to provide a mean estimate of \({E} \left[\ell\left(\omega\right)\right]\) close to \(1\), the importance sampling estimates have smaller variance.
Next, we present distributions of estimates for \(\hat{E} \left[L\left(\omega^t\right)\right]\), in cases for \(T=1, 5, 10, 20\).
def simulate_multiple_T(key, p_a, p_b, q_a, q_b, T_values,
N_simu, N_samples=10000):
"""Run simulate once per T, returning a dict keyed by T."""
keys = jax.random.split(key, len(T_values))
return {T: simulate(keys[i], p_a, p_b, q_a, q_b,
N_simu, T, N_samples)
for i, T in enumerate(T_values)}
The next function draws the histograms that we use to compare the two methods.
def plot_estimates(T_values, mc, imp, imp_label, n_rows=1):
"""Compare Monte Carlo and importance sampling estimates."""
n_cols = len(T_values) // n_rows
fig, axs = plt.subplots(n_rows, n_cols,
figsize=(14, 5 * n_rows))
μ_range = np.linspace(0, 2, 100)
for ax, T in zip(np.ravel(axs), T_values):
μ_L_g, μ_L_h = np.asarray(mc[T]), np.asarray(imp[T])
ax.set_xlabel('$μ_L$')
ax.set_ylabel('frequency')
ax.set_title(f'$T$={T}')
ax.hist(μ_L_g, bins=μ_range,
color='r', alpha=0.5, label='$g$ generating')
ax.hist(μ_L_h, bins=μ_range,
color='b', alpha=0.5, label=imp_label)
ax.legend(loc=4)
# Summarize each distribution in an upper corner
for μ_L, color, x, ha in ((μ_L_g, 'r', 0.02, 'left'),
(μ_L_h, 'b', 0.98, 'right')):
ax.text(x, 0.98, transform=ax.transAxes, ha=ha,
va='top', color=color, fontsize=9,
s=r'$\hat{μ}$=' + f'{np.mean(μ_L):.3g}' +
'\n' + 'med=' + f'{np.median(μ_L):.3g}' +
'\n' + r'$\hat{σ}^2$=' + f'{np.var(μ_L):.3g}')
plt.show()
The Monte Carlo estimates do not depend on the importance distribution, so we compute them once and reuse them in each of the figures below.
T_values = [1, 5, 10, 20]
mc = simulate_multiple_T(jax.random.key(6), g_a, g_b,
g_a, g_b, T_values, N_simu)
imp_h1 = simulate_multiple_T(jax.random.key(7), g_a, g_b,
h_a, h_b, T_values, N_simu)
plot_estimates(T_values, mc, imp_h1, '$h$ generating',
n_rows=2)
Fig. 25.4 Monte Carlo and importance sampling estimates#
The simulation exercises above show that the importance sampling estimates stay centered on \(1\) for every \(T\), while the distribution of the standard Monte Carlo estimates drifts steadily to the left as \(T\) increases.
It is worth being careful about what is going wrong here.
The Monte Carlo estimator is unbiased in population — its expectation is exactly \(1\) for every \(T\).
What deteriorates as \(T\) grows is the shape of its sampling distribution: almost all of the mass collapses towards zero, while the expectation is sustained by rare, very large outliers.
This is why each panel reports a median alongside the mean.
The median falls steadily in \(T\), whereas the reported mean \(\hat{μ}\) is an unreliable statistic — recall that the underlying variance is infinite, so the \(\hat{σ}^2\) printed in each panel estimates a quantity that does not exist, and a different seed can move \(\hat{μ}\) a long way.
Importance sampling with \(h\) removes the infinite variance, but it does not stop the variance from growing with \(T\).
Drawing a path from \(h\) and reweighting gives the path the value \(\prod_{i=1}^T \ell(\omega_i) g(\omega_i)/h(\omega_i) = \prod_{i=1}^T f(\omega_i)/h(\omega_i)\).
Its second moment is \(\left(\int_0^1 f^2/h \, d\omega\right)^T = \left(1 + \chi^2(f\|h)\right)^T\), where \(\chi^2(f\|h) = \int_0^1 (f-h)^2/h \, d\omega\) is the \(\chi^2\) divergence of \(f\) from \(h\) (see Statistical Divergence Measures).
For \(h = \) Beta\((0.5, 0.5)\), each weight \(f(\omega)/h(\omega) = \pi \sqrt{\omega(1-\omega)}\) is at most \(\pi/2\), and \(\int_0^1 f^2/h \, d\omega = \pi^2/8 \approx 1.234\).
The variance of a single path’s contribution, \((\pi^2/8)^T - 1\), is finite but grows geometrically: it equals \(66\) at \(T = 20\) and about \(3.6 \times 10^4\) at \(T = 50\).
25.6. Choosing a sampling distribution#
Above, we arbitrarily chose \(h = Beta(0.5,0.5)\) as the importance distribution.
Is there an optimal importance distribution?
An importance sampling estimator has zero variance when the importance density is proportional to the integrand, here \(\ell g\).
Since \(\ell g = f\), the zero-variance importance density for a single draw is \(f\) itself; for a path of length \(T\) it is the product density \(f(\omega_1) f(\omega_2) \cdots f(\omega_T)\).
To use it we need only be able to draw from \(f\); we do not need to know in advance that \(E \left[ L\left(\omega^t\right) \right] = 1\).
Thus, suppose that we simply use \(h = f\).
When estimating the mean of the likelihood ratio (T=1), we get:
μ_L_f = simulate(jax.random.key(8), g_a, g_b,
params.F_a, params.F_b, N_simu)
# importance sampling (mean and variance)
jnp.mean(μ_L_f), jnp.var(μ_L_f)
(Array(1., dtype=float64), Array(0., dtype=float64))
We could also use other distributions as our importance distribution.
Below we choose just a few and compare their sampling properties.
a_list = [0.5, 1., 2.]
b_list = [0.5, 1.2, 5.]
w_range = np.linspace(1e-5, 1-1e-5, 1000)
plt.plot(w_range, g(w_range),
lw=2, label=f'g=Beta({g_a}, {g_b})')
plt.plot(w_range, beta.pdf(w_range, a_list[0], b_list[0]),
lw=2, label=f'$h_1$=Beta({a_list[0]},{b_list[0]})')
plt.plot(w_range, beta.pdf(w_range, a_list[1], b_list[1]),
lw=2, label=f'$h_2$=Beta({a_list[1]},{b_list[1]})')
plt.plot(w_range, beta.pdf(w_range, a_list[2], b_list[2]),
lw=2, label=f'$h_3$=Beta({a_list[2]},{b_list[2]})')
plt.legend()
plt.ylim([0., 3.])
plt.show()
Fig. 25.5 Comparison of importance sampling distributions#
We consider two additional distributions.
As a reminder \(h_1\) is the original \(Beta(0.5,0.5)\) distribution that we used above.
\(h_2\) is the \(Beta(1,1.2)\) distribution.
Note how \(h_2\) has a similar shape to \(g\) at higher values of \(\omega\) but more mass at lower values.
Our hunch is that \(h_2\) should be a good importance sampling distribution.
\(h_3\) is the \(Beta(2,5)\) distribution.
Note how \(h_3\) has almost no mass at values very close to 0 and at values close to 1.
Our hunch is that \(h_3\) will be a poor importance sampling distribution.
The variance calculation above lets us make these hunches precise.
Drawing from \(h\) and reweighting gives each observation the value \(\ell\left(\omega\right) g\left(\omega\right) / h\left(\omega\right) = f\left(\omega\right) / h\left(\omega\right)\), so the estimator has finite variance exactly when
For a Beta\((a, b)\) importance density \(h\) and the uniform \(f\),
when \(a < 2\) and \(b < 2\), and the integral is infinite otherwise.
The integral equals \(\pi^2/8 \approx 1.234\) for \(h_1\), \(1/0.96 \approx 1.042\) for \(h_2\), and \(1\) for \(h = f\).
The proposal \(h_2\) is good because it is close to \(f\).
For \(h_3\) it diverges at both endpoints, since \(h_3\) vanishes like \(\omega\) at the origin and like \(\left(1-\omega\right)^4\) at one.
So importance sampling with \(h_3\) has infinite variance, just as standard Monte Carlo does — which is what makes it a poor choice.
We first simulate and plot the distribution of estimates for \(\hat{E} \left[L\left(\omega^t\right)\right]\) using \(h_2\) as the importance sampling distribution.
T_values_h2 = [1, 20]
imp_h2 = simulate_multiple_T(jax.random.key(9), g_a, g_b,
a_list[1], b_list[1],
T_values_h2, N_simu)
plot_estimates(T_values_h2, mc, imp_h2, '$h_2$ generating')
Fig. 25.6 Estimates using importance distribution \(h_2\)#
Our simulations suggest that indeed \(h_2\) is a quite good importance sampling distribution for our problem.
Even at \(T=20\), the mean is very close to \(1\) and the variance is small.
T_values_h3 = [1, 20]
imp_h3 = simulate_multiple_T(jax.random.key(10), g_a, g_b,
a_list[2], b_list[2],
T_values_h3, N_simu)
plot_estimates(T_values_h3, mc, imp_h3, '$h_3$ generating')
Fig. 25.7 Estimates using importance distribution \(h_3\)#
However, \(h_3\) is evidently a poor importance sampling distribution for our problem, with a mean estimate far away from \(1\) for \(T = 20\).
Because \(h_3 > 0\) on \((0,1)\), the importance sampling estimator with \(h_3\) is unbiased for every \(T\), just as the Monte Carlo estimator under \(g\) is.
The trouble lies in the shape of its sampling distribution.
Even at \(T = 1\), the estimates under \(h_3\) are more dispersed than the Monte Carlo estimates under \(g\), and the median estimate lies farther below \(1\).
Near \(\omega = 1\), \(h_3\) vanishes like \((1-\omega)^4\) while \(f = 1\), so the weight \(f/h_3\) has tail index \(5/4\) under \(h_3\), smaller than the tail index \(3/2\) of \(f/g\) under \(g\).
The error of a sample mean of \(N\) draws therefore shrinks like \(N^{-1/5}\) under \(h_3\) and like \(N^{-1/3}\) under \(g\) (Exercise 25.2).
Thus, our simulations suggest that for our problem we would be better off simply using Monte Carlo approximations under \(g\) than using \(h_3\) as an importance sampling distribution.
25.7. Exercises#
Exercise 25.1
Consider estimating \(E^g\left[L\left(\omega^T\right)\right] = 1\) by drawing paths of length \(T\) from a Beta\((a, b)\) importance density \(q\).
Show that each path contributes \(\prod_{i=1}^T f(\omega_i)/q(\omega_i)\) to the estimator and that the second moment of this contribution is \(M(a,b)^T\), where \(M(a,b) = \int_0^1 f^2/q \, d\omega = B(a,b) B(2-a, 2-b)\) when \(a, b < 2\) and \(M(a,b) = +\infty\) otherwise.
Use the Cauchy–Schwarz inequality to show that \(M(a,b) \geq 1\), with equality only at \((a,b) = (1,1)\). Plot contours of \(\log M(a,b)\) on \((0,2) \times (0,2)\).
The effective sample size of \(N\) weighted paths is approximately \(N / M(a,b)^T\). Plot the ratio of the effective sample size to \(N\) against \(T\) for \(h_1\) and \(h_2\). How many paths does \(h_1\) need for an effective sample size of \(1000\) at \(T = 20\) and at \(T = 50\)?
Solution
Part 1. Each path contributes \(\prod_{i=1}^T \ell(\omega_i) g(\omega_i)/q(\omega_i) = \prod_{i=1}^T f(\omega_i)/q(\omega_i)\).
Because the \(\omega_i\)’s are IID draws from \(q\), the second moment is \(\left(E^q\left[(f/q)^2\right]\right)^T = M(a,b)^T\) with \(M(a,b) = \int_0^1 f^2/q \, d\omega\).
Since \(f = 1\), \(M(a,b) = B(a,b) \int_0^1 \omega^{1-a} (1-\omega)^{1-b} d\omega\), which equals \(B(a,b) B(2-a, 2-b)\) when \(a < 2\) and \(b < 2\) and is infinite otherwise.
Part 2. By the Cauchy–Schwarz inequality,
with equality if and only if \(f/\sqrt{q}\) is proportional to \(\sqrt{q}\), that is, \(q = f\).
The variance of a single weight is \(M(a,b) - 1 = \chi^2(f\|q)\).
from scipy.special import betaln
def M(a, b):
"""Second moment of the importance weight f/q for q = Beta(a, b)."""
if a < 2 and b < 2:
return np.exp(betaln(a, b) + betaln(2 - a, 2 - b))
return np.inf
grid = np.linspace(0.05, 1.95, 191)
log_M = np.array([[np.log(M(a, b)) for a in grid] for b in grid])
fig, ax = plt.subplots(figsize=(7, 6))
cs = ax.contourf(grid, grid, np.minimum(log_M, 2),
levels=np.linspace(0, 2, 21), cmap='viridis')
fig.colorbar(cs, label=r'$\log M(a, b)$ (capped at 2)')
lines = ax.contour(grid, grid, log_M, levels=[0.01, 0.05, 0.2, 0.5, 1],
colors='w', linewidths=0.8)
ax.clabel(lines, fmt='%.2f', fontsize=8)
for (a, b), name in zip([(0.5, 0.5), (1.0, 1.2), (1.0, 1.0)],
['$h_1$', '$h_2$', '$f$']):
ax.plot(a, b, 'wo')
ax.annotate(name, (a, b), xytext=(5, 5), textcoords='offset points',
color='w')
ax.set_xlabel('$a$')
ax.set_ylabel('$b$')
plt.show()
i_b, i_a = np.unravel_index(np.argmin(log_M), log_M.shape)
print(f"Minimum of M on the grid at (a, b) = ({grid[i_a]:.2f}, {grid[i_b]:.2f})")
print(f"M(h_1) = {M(0.5, 0.5):.4f}, M(h_2) = {M(1, 1.2):.4f}")
Minimum of M on the grid at (a, b) = (1.00, 1.00)
M(h_1) = 1.2337, M(h_2) = 1.0417
Part 3.
T_grid = np.arange(1, 51)
plt.figure(figsize=(8, 5))
for (a, b), name in zip([(0.5, 0.5), (1.0, 1.2)], ['$h_1$', '$h_2$']):
plt.semilogy(T_grid, M(a, b)**(-T_grid), lw=2, label=name)
plt.xlabel('$T$')
plt.ylabel('effective sample size / $N$')
plt.legend()
plt.show()
for T in [20, 50]:
print(f"T = {T}: h_1 needs N = {1000 * M(0.5, 0.5)**T:,.0f} paths, "
f"h_2 needs N = {1000 * M(1, 1.2)**T:,.0f} paths")
T = 20: h_1 needs N = 66,711 paths, h_2 needs N = 2,262 paths
T = 50: h_1 needs N = 36,348,619 paths, h_2 needs N = 7,699 paths
The effective sample size under \(h_1\) shrinks by the factor \(8/\pi^2\) with each additional period.
Exercise 25.2
Let \(\omega \sim g\). Show that \(\Pr\{\ell(\omega) > x\}\) is approximately proportional to \(x^{-3/2}\) for large \(x\), so that \(\ell(\omega)\) has tail index \(3/2\) under \(g\).
Let \(\omega \sim h_3 = \) Beta\((2, 5)\). Show that the weight \(f(\omega)/h_3(\omega)\) has tail index \(5/4\) under \(h_3\). Which endpoint of \((0,1)\) determines the tail index?
Confirm both tail indices with the Hill estimator [Hill, 1975], which uses the \(k\) largest of \(n\) draws \(X_{(1)} \geq X_{(2)} \geq \cdots \geq X_{(k+1)}\) to estimate the tail index by
\[ \hat \alpha_k = \left(\frac{1}{k} \sum_{i=1}^k \log \frac{X_{(i)}}{X_{(k+1)}}\right)^{-1} . \]Set \(T = 1\). For \(N = 10^2, 10^3, 10^4, 10^5\), compute \(500\) replicate sample means of \(f/q\) for \(N\) draws from \(q\) when \(q\) is \(g\) (standard Monte Carlo), \(h_2\), and \(h_3\). Plot the interquartile range of the replicate means against \(N\) on log-log scales, and compare the slopes with \(-1/3\), \(-1/2\), and \(-1/5\).
Solution
Part 1. Near the origin, \(g(\omega) \approx k \omega^2\) with \(k = 1/B(3, 1.2)\).
For large \(x\), \(\ell(\omega) = 1/g(\omega) > x\) occurs mainly when \(\omega < (kx)^{-1/2}\).
Since \(\Pr_g\{\omega < u\} \approx k u^3/3\) for small \(u\),
Near \(\omega = 1\), \(g\) vanishes like \((1-\omega)^{0.2}\), which contributes a tail proportional to \(x^{-6}\), too thin to matter.
Part 2. The density \(h_3(\omega) = 30\, \omega (1-\omega)^4\) vanishes at both endpoints.
Near \(\omega = 0\), \(1/h_3 > x\) occurs when \(\omega < 1/(30x)\), which has probability about \(15/(30x)^2\), proportional to \(x^{-2}\).
Near \(\omega = 1\), \(1/h_3 > x\) occurs when \(1 - \omega < (30x)^{-1/4}\), which has probability about \(6 (30x)^{-5/4}\), proportional to \(x^{-5/4}\).
The heavier tail, with index \(5/4\), comes from the right endpoint.
By the generalized central limit theorem, a tail index \(\alpha \in (1,2)\) makes the error of a sample mean shrink like \(N^{1/\alpha - 1}\): \(N^{-1/3}\) for \(\alpha = 3/2\) and \(N^{-1/5}\) for \(\alpha = 5/4\).
Part 3.
def log_beta_pdf(w, a, b):
"""Log density of Beta(a, b)."""
return (a - 1) * np.log(w) + (b - 1) * np.log1p(-w) - betaln(a, b)
def hill(x, k):
"""Hill estimator of the tail index from the k largest values of x."""
x_sorted = np.sort(x)[::-1]
return 1 / np.mean(np.log(x_sorted[:k] / x_sorted[k]))
rng = np.random.default_rng(0)
n = 1_000_000
w_g = rng.beta(g_a, g_b, n)
w_h3 = rng.beta(2, 5, n)
weights = {'f/g under g': np.exp(-log_beta_pdf(w_g, g_a, g_b)),
'f/h_3 under h_3': np.exp(-log_beta_pdf(w_h3, 2, 5))}
for name, x in weights.items():
estimates = [f"{hill(x, k):.3f}" for k in (500, 1000, 2000, 5000)]
print(f"{name}: Hill estimates for k = 500, 1000, 2000, 5000: {estimates}")
f/g under g: Hill estimates for k = 500, 1000, 2000, 5000: ['1.483', '1.476', '1.497', '1.514']
f/h_3 under h_3: Hill estimates for k = 500, 1000, 2000, 5000: ['1.248', '1.265', '1.267', '1.271']
The Hill estimates are close to \(3/2\) and \(5/4\).
Part 4.
N_values = np.array([10**2, 10**3, 10**4, 10**5])
R = 500
proposals = {'$g$ (standard Monte Carlo)': (g_a, g_b, 1/3),
'$h_2$': (1.0, 1.2, 1/2),
'$h_3$': (2.0, 5.0, 1/5)}
plt.figure(figsize=(8, 5))
for name, (a, b, rate) in proposals.items():
iqr = []
for N in N_values:
means = np.empty(R)
for r in range(R):
w = np.clip(rng.beta(a, b, N), 1e-16, 1 - 1e-16)
means[r] = np.mean(np.exp(-log_beta_pdf(w, a, b)))
q75, q25 = np.quantile(means, [0.75, 0.25])
iqr.append(q75 - q25)
slope = np.polyfit(np.log(N_values), np.log(iqr), 1)[0]
line, = plt.loglog(N_values, iqr, 'o-',
label=f'{name}: slope {slope:.2f}')
plt.loglog(N_values, iqr[0] * (N_values / N_values[0])**(-rate),
'--', color=line.get_color(), lw=1)
plt.xlabel('$N$')
plt.ylabel('interquartile range of replicate means')
plt.legend()
plt.show()
The dashed lines have slopes \(-1/3\), \(-1/2\), and \(-1/5\).
The fitted slopes are close to them.
Exercise 25.3
In this section of Likelihood Ratio Processes, the probability \(\alpha_T = \Pr_f\{L(\omega^T) < 1\}\) of a type I error decays exponentially in \(T\), so standard Monte Carlo cannot estimate it for large \(T\).
Let \(\phi^*\) minimize the Chernoff coefficient \(\rho(\phi) = \int_0^1 f^\phi g^{1-\phi} d\omega\), and let \(q = f^{\phi^*} g^{1-\phi^*}/\rho(\phi^*)\) be the tilted density.
Show that \(q\) is the Beta\((\phi^* F_a + (1-\phi^*) G_a, \ \phi^* F_b + (1-\phi^*) G_b)\) density, and compute \(\phi^*\) and the parameters of \(q\).
Show that
\[ \alpha_T = \rho(\phi^*)^T \, E^q\left[\mathbf{1}\{L(\omega^T) < 1\} \, L(\omega^T)^{1-\phi^*}\right] . \]For \(T = 30, 60, 100\), estimate \(\alpha_T\) by standard Monte Carlo with \(N = 20{,}000\) paths drawn from \(f\) and by importance sampling with \(N = 20{,}000\) paths drawn from \(q\). Report the relative standard error of the importance sampling estimate.
Solution
Part 1. For Beta densities,
which is proportional to the stated Beta density.
Part 2. Since \(f/q = \rho(\phi^*) (f/g)^{1-\phi^*}\),
On the event \(L(\omega^T) < 1\), the weight \(\rho(\phi^*)^T L(\omega^T)^{1-\phi^*}\) is at most \(\rho(\phi^*)^T = e^{-C(f,g) T}\), so the estimator has a small variance.
At \(\phi^*\), the derivative of \(\log \rho(\phi)\), which equals \(E^q \log \ell(\omega)\), is zero.
Thus \(\log L(\omega^T)\) has mean zero under \(q\), and the rare event \(L(\omega^T) < 1\) under \(f\) occurs about half the time under \(q\).
Choosing \(q\) in this way is called exponential tilting [Bucklew, 2004].
Part 3.
from scipy.optimize import minimize_scalar
F_a, F_b = params.F_a, params.F_b
def log_ρ(ϕ):
"""Log Chernoff coefficient for f = Beta(F_a, F_b) and g = Beta(g_a, g_b)."""
return (betaln(ϕ * F_a + (1 - ϕ) * g_a, ϕ * F_b + (1 - ϕ) * g_b)
- ϕ * betaln(F_a, F_b) - (1 - ϕ) * betaln(g_a, g_b))
res = minimize_scalar(log_ρ, bounds=(0, 1), method='bounded')
ϕ_star, log_ρ_star = res.x, res.fun
q_a = ϕ_star * F_a + (1 - ϕ_star) * g_a
q_b = ϕ_star * F_b + (1 - ϕ_star) * g_b
print(f"ϕ* = {ϕ_star:.4f}, C(f, g) = {-log_ρ_star:.4f}, "
f"q = Beta({q_a:.4f}, {q_b:.4f})")
rng = np.random.default_rng(1)
N = 20_000
for T in [30, 60, 100]:
# Standard Monte Carlo under f
w = rng.uniform(size=(N, T))
log_L = np.sum(log_beta_pdf(w, F_a, F_b) - log_beta_pdf(w, g_a, g_b), axis=1)
α_mc = np.mean(log_L < 0)
# Importance sampling under the tilted density q
w = np.clip(rng.beta(q_a, q_b, (N, T)), 1e-300, 1 - 1e-16)
log_L = np.sum(log_beta_pdf(w, F_a, F_b) - log_beta_pdf(w, g_a, g_b), axis=1)
X = np.where(log_L < 0, np.exp(T * log_ρ_star + (1 - ϕ_star) * log_L), 0.0)
α_is = np.mean(X)
rel_se = np.std(X) / np.sqrt(N) / α_is
print(f"T = {T:>3}: standard MC = {α_mc:.2e}, "
f"importance sampling = {α_is:.3e} (relative s.e. {rel_se:.3f})")
ϕ* = 0.5969, C(f, g) = 0.1212, q = Beta(1.8063, 1.0806)
T = 30: standard MC = 4.15e-03, importance sampling = 4.401e-03 (relative s.e. 0.011)
T = 60: standard MC = 1.00e-04, importance sampling = 8.448e-05 (relative s.e. 0.013)
T = 100: standard MC = 0.00e+00, importance sampling = 5.143e-07 (relative s.e. 0.015)
At \(T = 30\) the two estimates agree.
At \(T = 60\), standard Monte Carlo rests on a handful of paths with \(L(\omega^T) < 1\), while importance sampling has a relative standard error near \(1\) percent.
At \(T = 100\), standard Monte Carlo records no paths with \(L(\omega^T) < 1\), while importance sampling still delivers an estimate with a relative standard error near \(1.5\) percent.