22. Likelihood Ratio Processes#
In addition to what’s in Anaconda, this lecture will need the following libraries:
!pip install --upgrade quantecon
22.1. Overview#
This lecture describes likelihood ratio processes and some of their uses.
We’ll study the same setting that is also used in this lecture on exchangeability.
Among the things that we’ll learn are
How a likelihood ratio process is a key ingredient in frequentist hypothesis testing
How a receiver operating characteristic curve summarizes information about a false alarm probability and power in frequentist hypothesis testing
How a statistician can combine frequentist probabilities of type I and type II errors to form posterior probabilities of mistakes in a model selection or in an individual-classification problem
How to use a Kullback-Leibler divergence to quantify the difference between two probability distributions with the same support
How during World War II the United States Navy devised a decision rule for doing quality control on lots of ammunition, a topic that sets the stage for this lecture
A peculiar property of likelihood ratio processes
Let’s start by importing some Python tools.
import matplotlib.pyplot as plt
import numpy as np
from numba import vectorize, jit
from math import gamma
from scipy.integrate import quad
from scipy.optimize import brentq, minimize_scalar
from scipy.stats import beta as beta_dist
from scipy.special import betaln, digamma
import pandas as pd
from IPython.display import display, Math
import quantecon as qe
rng = np.random.default_rng()
22.2. Likelihood Ratio Process#
A nonnegative random variable \(W\) has one of two probability density functions, either \(f\) or \(g\).
Before the beginning of time, nature once and for all decides whether she will draw a sequence of IID draws from either \(f\) or \(g\).
We let \(q\) be the density that nature chose once and for all, so that \(q\) is either \(f\) or \(g\), permanently.
Nature knows which density it permanently draws from, but we the observers do not.
We know both \(f\) and \(g\) but we don’t know which density nature chose.
But we want to know.
To do that, we use observations.
We observe a sequence \(\{w_t\}_{t=1}^T\) of \(T\) IID draws that we know came from either \(f\) or \(g\).
We want to use these observations to infer whether nature chose \(f\) or \(g\).
A likelihood ratio process is a useful tool for this task.
To begin, we define a key component of a likelihood ratio process, namely, the time \(t\) likelihood ratio as the random variable
We assume that \(f\) and \(g\) both put positive probabilities on the same intervals of possible realizations of the random variable \(W\).
That means that under the \(g\) density, \(\ell (w_t)= \frac{f\left(w_{t}\right)}{g\left(w_{t}\right)}\) is a nonnegative random variable with mean \(1\).
A likelihood ratio process for sequence \(\left\{ w_{t}\right\} _{t=1}^{\infty}\) is defined as
where \(w^t=\{ w_1,\dots,w_t\}\) is a history of observations up to and including time \(t\).
Sometimes for shorthand we’ll write \(L_t = L(w^t)\).
Notice that the likelihood process satisfies the recursion
The likelihood ratio and its logarithm are key tools for making inferences using a classic frequentist approach due to Neyman and Pearson [Neyman and Pearson, 1933].
To help us appreciate how things work, the following Python code evaluates \(f\) and \(g\) as two different Beta distributions, then computes and simulates an associated likelihood ratio process by generating a sequence \(w^t\) from one of the two probability distributions, for example, a sequence of IID draws from \(g\).
# Parameters for the two Beta distributions
F_a, F_b = 1, 1
G_a, G_b = 3, 1.2
@vectorize
def p(x, a, b):
"""Beta distribution density function."""
r = gamma(a + b) / (gamma(a) * gamma(b))
return r * x** (a-1) * (1 - x) ** (b-1)
f = jit(lambda x: p(x, F_a, F_b))
g = jit(lambda x: p(x, G_a, G_b))
def create_beta_density(a, b):
"""Create a beta density function with specified parameters."""
return jit(lambda x: p(x, a, b))
def likelihood_ratio(w, f_func, g_func):
"""Compute likelihood ratio for observation(s) w."""
return f_func(w) / g_func(w)
@jit
def simulate_likelihood_ratios(a, b, f_func, g_func, rng, T=50, N=500):
"""
Generate N sets of T observations of the likelihood ratio.
"""
l_arr = np.empty((N, T))
for i in range(N):
for j in range(T):
w = rng.beta(a, b)
l_arr[i, j] = f_func(w) / g_func(w)
return l_arr
def simulate_sequences(distribution, f_func, g_func, rng,
F_params=(1, 1), G_params=(3, 1.2), T=50, N=500):
"""
Generate N sequences of T observations from specified distribution.
"""
if distribution == 'f':
a, b = F_params
elif distribution == 'g':
a, b = G_params
else:
raise ValueError("distribution must be 'f' or 'g'")
l_arr = simulate_likelihood_ratios(a, b, f_func, g_func, rng, T, N)
l_seq = np.cumprod(l_arr, axis=1)
return l_arr, l_seq
def plot_likelihood_paths(l_seq, title="Likelihood ratio paths",
ylim=None, n_paths=None):
"""Plot likelihood ratio paths."""
N, T = l_seq.shape
n_show = n_paths or min(N, 100)
plt.figure(figsize=(10, 6))
for i in range(n_show):
plt.plot(np.arange(1, T+1), l_seq[i, :], color='b', lw=0.8, alpha=0.5)
if ylim:
plt.ylim(ylim)
plt.title(title)
plt.xlabel('t')
plt.ylabel('$L(w^t)$')
plt.show()
22.3. Nature permanently draws from density g#
We first simulate the likelihood ratio process when nature permanently draws from \(g\).
# Simulate when nature draws from g
l_arr_g, l_seq_g = simulate_sequences('g', f, g, rng, (F_a, F_b), (G_a, G_b))
plot_likelihood_paths(l_seq_g,
title="$L(w^{t})$ paths when nature draws from g",
ylim=[0, 3])
Evidently, as sample length \(T\) grows, most probability mass shifts toward zero.
To see this more clearly, we plot over time the fraction of paths \(L\left(w^{t}\right)\) that fall in the interval \(\left[0, 0.01\right]\).
N, T = l_arr_g.shape
plt.plot(np.arange(1, T+1), np.sum(l_seq_g <= 0.01, axis=0) / N)
plt.xlabel('t')
plt.show()
Despite the evident convergence of most probability mass to a very small interval near \(0\), the unconditional mean of \(L\left(w^t\right)\) under probability density \(g\) is identically \(1\) for all \(t\).
To verify this assertion, first notice that as mentioned earlier the unconditional mean \(E\left[\ell \left(w_{t}\right)\bigm|q=g\right]\) is \(1\) for all \(t\):
which immediately implies
Because \(L(w^t) = \ell(w_t) L(w^{t-1})\) and \(\{w_t\}_{t=1}^\infty\) is an IID sequence, we have
for any \(t \geq 1\).
Mathematical induction implies \(E\left[L\left(w^{t}\right)\bigm|q=g\right]=1\) for all \(t \geq 1\).
22.4. Peculiar property#
How can \(E\left[L\left(w^{t}\right)\bigm|q=g\right]=1\) possibly be true when most probability mass of the likelihood ratio process is piling up near \(0\) as \(t \rightarrow + \infty\)?
The answer is that as \(t \rightarrow + \infty\), the distribution of \(L_t\) becomes more and more fat-tailed: enough mass shifts to larger and larger values of \(L_t\) to make the mean of \(L_t\) continue to be one despite most of the probability mass piling up near \(0\).
Under \(g\), the process \(\{L_t\}\) is a nonnegative martingale with mean \(1\).
By the strong law of large numbers, \(\frac{1}{t} \log L_t \rightarrow E_g \log \ell(w) = -D_{KL}(g\|f) < 0\) almost surely, so \(L_t \rightarrow 0\) almost surely.
Thus \(E[L_t \mid q=g] = 1\) for every \(t\) while \(E[\lim_{t \rightarrow \infty} L_t \mid q=g] = 0\).
The martingale \(\{L_t\}\) converges almost surely but not in mean, because it is not uniformly integrable.
Merging of Opinions: The Blackwell–Dubins Theorem studies this likelihood-ratio martingale.
Kakutani’s theorem there states that when \(\int \sqrt{f g} \, dw < 1\), the distributions of infinite sequences of IID draws from \(f\) and from \(g\) are mutually singular, which is another way of saying that \(L_t \rightarrow 0\) almost surely under \(g\).
To illustrate this peculiar property, we simulate many paths and calculate the unconditional mean of \(L\left(w^t\right)\) by averaging across these many paths at each \(t\).
l_arr_g, l_seq_g = simulate_sequences('g',
f, g, rng, (F_a, F_b), (G_a, G_b), N=50000)
It would be useful to use simulations to verify that unconditional means \(E\left[L\left(w^{t}\right)\right]\) equal unity by averaging across sample paths.
But it would be too computer-time-consuming for us to do that here simply by applying a standard Monte Carlo simulation approach.
The reason is that the distribution of \(L\left(w^{t}\right)\) is extremely skewed for large values of \(t\).
Because the probability density in the right tail is close to \(0\), it just takes too much computer time to sample enough points from the right tail.
We explain the problem in more detail in this lecture.
There we describe an alternative way to compute the mean of a likelihood ratio by computing the mean of a different random variable by sampling from a different probability distribution.
22.5. Nature permanently draws from density f#
Now suppose that before time \(0\) nature permanently decided to draw repeatedly from density \(f\).
While the mean of the likelihood ratio \(\ell \left(w_{t}\right)\) under density \(g\) is \(1\), its mean under the density \(f\) exceeds one.
To see this, we compute
For the \(f\) and \(g\) in our example, \(E\left[\ell \left(w_{t}\right)\bigm|q=f\right]\) is infinite.
Because \(f\) is the uniform density, \(E\left[\ell \left(w_{t}\right)\bigm|q=f\right] = \int_0^1 1/g(w) \, dw\).
As \(w \rightarrow 0\), \(g(w)\) vanishes like \(w^{G_a - 1} = w^2\), so the integral diverges.
Exercise 22.4 shows that \(E\left[\ell \left(w_{t}\right)\bigm|q=f\right] = B(G_a, G_b) B(2-G_a, 2-G_b)\) when \(G_a < 2\) and \(G_b < 2\), where \(B\) is the Beta function, and that it is infinite otherwise.
Because the \(w_t\)’s are IID, \(E\left[L\left(w^{t}\right)\bigm|q=f\right] = +\infty\) for every \(t \geq 1\).
A sample average of \(L\left(w^{t}\right)\) across simulated paths therefore does not settle down as the number of paths grows.
Mean of a Likelihood Ratio Process encounters the same divergent integral as the infinite variance of \(\ell(w_t)\) under \(g\), since \(E\left[\ell \left(w_{t}\right)^2\bigm|q=g\right] = E\left[\ell \left(w_{t}\right)\bigm|q=f\right]\).
The logarithm of \(L\left(w^{t}\right)\) behaves well.
Under \(f\), \(E\left[\log L\left(w^{t}\right)\bigm|q=f\right] = t \, D_{KL}(f\|g)\), and by the strong law of large numbers \(\frac{1}{t} \log L\left(w^{t}\right) \rightarrow D_{KL}(f\|g) > 0\) almost surely, so \(L\left(w^{t}\right) \rightarrow + \infty\) almost surely.
# Simulate when nature draws from f
l_arr_f, l_seq_f = simulate_sequences('f', f, g, rng,
(F_a, F_b), (G_a, G_b), N=50000)
The next figure plots the mean and the median across paths of \(\log L\left(w^{t}\right)\), together with the line \(t \, D_{KL}(f\|g)\).
N, T = l_arr_f.shape
t_grid = np.arange(1, T+1)
KL_fg, _ = quad(lambda w: f(w) * np.log(f(w) / g(w)), 0, 1)
log_l_seq_f = np.log(l_seq_f)
plt.plot(t_grid, np.mean(log_l_seq_f, axis=0), label=r'mean of $\log L_t$')
plt.plot(t_grid, np.median(log_l_seq_f, axis=0), label=r'median of $\log L_t$')
plt.plot(t_grid, KL_fg * t_grid, 'k--', label=r'$t \, D_{KL}(f\|g)$')
plt.xlabel('t')
plt.legend()
plt.show()
print(f"D_KL(f||g) = {KL_fg:.4f}")
D_KL(f||g) = 0.7592
We also plot the probability that \(L\left(w^t\right)\) falls into the interval \([10000, \infty)\) as a function of time and watch how fast probability mass diverges to \(+\infty\).
plt.plot(t_grid, np.sum(l_seq_f > 10000, axis=0) / N)
plt.xlabel('t')
plt.show()
22.6. Likelihood ratio test#
We now describe how to employ the machinery of Neyman and Pearson [Neyman and Pearson, 1933] to test the hypothesis that history \(w^t\) is generated by repeated IID draws from density \(f\).
Denote \(q\) as the data generating process, so that \(q=f \text{ or } g\).
Upon observing a sample \(\{W_i\}_{i=1}^t\), we want to decide whether nature is drawing from \(g\) or from \(f\) by performing a (frequentist) hypothesis test.
We specify
Null hypothesis \(H_0\): \(q=f\),
Alternative hypothesis \(H_1\): \(q=g\).
Neyman and Pearson proved that the best way to test this hypothesis is to use a likelihood ratio test that takes the form:
accept \(H_0\) if \(L(W^t) > c\),
reject \(H_0\) if \(L(W^t) < c\),
where \(c\) is a given discrimination threshold.
Setting \(c =1\) is a common choice.
When \(f\) and \(g\) have equal prior probabilities, \(c = 1\) selects the model with the larger posterior probability, which minimizes the probability of selecting the wrong model.
We’ll discuss consequences of other choices of \(c\) below.
This test is best in the sense that it is most powerful.
To understand what this means, we have to define probabilities of two important events that allow us to characterize a test associated with a given threshold \(c\).
The two probabilities are:
Probability of a Type I error in which we reject \(H_0\) when it is true:
\[ \alpha \equiv \Pr\left\{ L\left(w^{t}\right)<c\mid q=f\right\} \]Probability of a Type II error in which we accept \(H_0\) when it is false:
\[ \beta \equiv \Pr\left\{ L\left(w^{t}\right)>c\mid q=g\right\} \]
The probability \(\alpha\) is also called the probability of false alarm or the significance level.
The probability \(1-\beta = \Pr\left\{ L\left(w^{t}\right)<c\mid q=g\right\}\) is called the probability of detection or the power of the test.
The Neyman-Pearson Lemma states that among all possible tests, a likelihood ratio test maximizes the probability of detection for a given probability of false alarm.
Another way to say the same thing is that among all possible tests, a likelihood ratio test maximizes power for a given significance level.
We want a small probability of false alarm and a large probability of detection.
With sample size \(t\) fixed, we can change our two probabilities by adjusting \(c\).
A troublesome “that’s life” fact is that these two probabilities move in the same direction as we vary the critical value \(c\).
Without specifying quantitative losses from making Type I and Type II errors, there is little that we can say about how we should trade off probabilities of the two types of mistakes.
We do know that increasing sample size \(t\) improves statistical inference.
The experiment that reveals \(w^t\) is more informative in Blackwell’s sense than the experiment that reveals \(w^s\) for \(t > s\), because a statistician can always discard the last \(t-s\) observations.
When there are only two hypotheses, one experiment is more informative than another in Blackwell’s sense if and only if its receiver operating characteristic curve, defined below, lies everywhere above the other’s.
Blackwell’s Theorem on Comparing Experiments describes Blackwell’s ordering of experiments; its section on Bayesian likelihood-ratio learning applies it to the setting of this lecture.
Below we plot some informative figures that illustrate this.
We also present a classical frequentist method for choosing a sample size \(t\).
Let’s start with a case in which we fix the threshold \(c\) at \(1\).
c = 1
Below we plot empirical distributions of logarithms of the cumulative likelihood ratios simulated above, which are generated by either \(f\) or \(g\).
Taking logarithms has no effect on calculating the probabilities because the log is a monotonic transformation.
As \(t\) increases, the probabilities of making Type I and Type II errors both decrease, which is good.
This is because most of the probability mass of log\((L(w^t))\) moves toward \(-\infty\) when \(g\) is the data generating process, while log\((L(w^t))\) goes to \(\infty\) when data are generated by \(f\).
That disparate behavior of log\((L(w^t))\) under \(f\) and \(g\) is what makes it possible eventually to distinguish \(q=f\) from \(q=g\).
def plot_log_histograms(l_seq_f, l_seq_g, c=1, time_points=[1, 7, 14, 21]):
"""Plot log likelihood ratio histograms."""
fig, axs = plt.subplots(2, 2, figsize=(12, 8))
for i, t in enumerate(time_points):
nr, nc = i // 2, i % 2
axs[nr, nc].axvline(np.log(c), color="k", ls="--")
# Column t-1 of l_seq holds L(w^t)
hist_f, x_f = np.histogram(np.log(l_seq_f[:, t-1]), 200, density=True)
hist_g, x_g = np.histogram(np.log(l_seq_g[:, t-1]), 200, density=True)
axs[nr, nc].plot(x_f[1:], hist_f, color="C0", label="dist under f")
axs[nr, nc].plot(x_g[1:], hist_g, color="C1", label="dist under g")
# Fill error regions
for j, (x, hist, label) in enumerate(
zip([x_f, x_g], [hist_f, hist_g],
["Type I error", "Type II error"])):
ind = x[1:] <= np.log(c) if j == 0 else x[1:] > np.log(c)
axs[nr, nc].fill_between(x[1:][ind], hist[ind], color=f"C{j}",
alpha=0.5, label=label)
axs[nr, nc].legend()
axs[nr, nc].set_title(f"t={t}")
plt.show()
plot_log_histograms(l_seq_f, l_seq_g, c=c)
In the above graphs,
the blue areas equal the probabilities \(\alpha\) of a type I error, up to the binning of the histograms
the orange areas equal the probabilities \(\beta\) of a type II error, up to the binning of the histograms
Because the logarithm is monotone, the event \(\{L_t < 1\}\) is the event \(\{\log L_t < 0\}\), so the area under the density of \(\log L_t\) to the left of \(0\) is \(\Pr\{L_t < 1\}\).
When we hold \(c\) fixed at \(c=1\), the following graph shows that
the probability of detection monotonically increases with increases in \(t\)
the probability of a false alarm monotonically decreases with increases in \(t\).
def compute_error_probabilities(l_seq_f, l_seq_g, c=1):
"""
Compute Type I and Type II error probabilities.
"""
N, T = l_seq_f.shape
# Type I error (false alarm) - reject H0 when true
PFA = np.array([np.sum(l_seq_f[:, t] < c) / N for t in range(T)])
# Type II error - accept H0 when false
beta = np.array([np.sum(l_seq_g[:, t] >= c) / N for t in range(T)])
# Probability of detection (power)
PD = np.array([np.sum(l_seq_g[:, t] < c) / N for t in range(T)])
return {
'alpha': PFA,
'beta': beta,
'PD': PD,
'PFA': PFA
}
def plot_error_probabilities(error_dict, T, c=1, title_suffix=""):
"""Plot error probabilities over time."""
t_grid = np.arange(1, T+1)
plt.figure(figsize=(10, 6))
plt.plot(t_grid, error_dict['PD'], label="Probability of detection")
plt.plot(t_grid, error_dict['PFA'], label="Probability of false alarm")
plt.xlabel("t")
plt.ylabel("Probability")
plt.title(f"Error Probabilities (c={c}){title_suffix}")
plt.legend()
plt.show()
error_probs = compute_error_probabilities(l_seq_f, l_seq_g, c=c)
N, T = l_seq_f.shape
plot_error_probabilities(error_probs, T, c)
For a given sample size \(t\), the threshold \(c\) uniquely pins down probabilities of both types of error.
If for a fixed \(t\) we now free up and move \(c\), we will sweep out the probability of detection as a function of the probability of false alarm.
This produces a receiver operating characteristic curve (ROC curve).
Below, we plot receiver operating characteristic curves for different sample sizes \(t\).
def plot_roc_curves(l_seq_f, l_seq_g, t_values=[1, 5, 9, 13], N=None):
"""Plot ROC curves for different sample sizes."""
if N is None:
N = l_seq_f.shape[0]
PFA = np.arange(0, 100, 1)
plt.figure(figsize=(10, 6))
for t in t_values:
# Column t-1 of l_seq holds L(w^t)
percentile = np.percentile(l_seq_f[:, t-1], PFA)
PD = [np.sum(l_seq_g[:, t-1] < p) / N for p in percentile]
plt.plot(PFA / 100, PD, label=f"t={t}")
plt.scatter(0, 1, label="perfect detection")
plt.plot([0, 1], [0, 1], color='k', ls='--', label="random detection")
plt.arrow(0.5, 0.5, -0.15, 0.15, head_width=0.03)
plt.text(0.35, 0.7, "better")
plt.xlabel("Probability of false alarm")
plt.ylabel("Probability of detection")
plt.legend()
plt.title("ROC Curve")
plt.show()
plot_roc_curves(l_seq_f, l_seq_g, t_values=range(1, 15, 4), N=N)
Notice that as \(t\) increases, we are assured a larger probability of detection and a smaller probability of false alarm associated with a given discrimination threshold \(c\).
For a given sample size \(t\), both \(\alpha\) and \(\beta\) change as we vary \(c\).
As we increase \(c\)
\(\alpha \equiv \Pr\left\{ L\left(w^{t}\right)<c\mid q=f\right\}\) increases
\(\beta \equiv \Pr\left\{ L\left(w^{t}\right)>c\mid q=g\right\}\) decreases
As \(t \rightarrow + \infty\), we approach the perfect detection curve that is indicated by a right angle hinging on the blue dot.
For a given sample size \(t\), the discrimination threshold \(c\) determines a point on the receiver operating characteristic curve.
It is up to the test designer to trade off probabilities of making the two types of errors.
But we know how to choose the smallest sample size to achieve given targets for the probabilities.
Typically, frequentists aim for a high probability of detection that respects an upper bound on the probability of false alarm.
Below we show an example in which we fix the probability of false alarm at \(0.05\).
The required sample size for making a decision is then determined by a target probability of detection, for example, \(0.9\), as depicted in the following graph.
PFA = 0.05
PD = np.empty(T)
for t in range(T):
c = np.percentile(l_seq_f[:, t], PFA * 100)
PD[t] = np.sum(l_seq_g[:, t] < c) / N
# PD[t] refers to sample size t+1
plt.plot(np.arange(1, T+1), PD)
plt.axhline(0.9, color="k", ls="--")
plt.xlabel("t")
plt.ylabel("Probability of detection")
plt.title(f"Probability of false alarm={PFA}")
plt.show()
The United States Navy evidently used a procedure like this to select a sample size \(t\) for doing quality control tests during World War II.
A Navy Captain who had been ordered to perform tests of this kind had doubts about it that he presented to Milton Friedman, as we describe in this lecture.
A Bayesian Formulation of Friedman and Wald’s Problem formulates the Captain’s problem as a Bayesian sequential decision problem, and Bayesian versus Frequentist Decision Rules compares the Navy’s frequentist rule with the Bayesian decision rule.
22.6.1. A third distribution \(h\)#
Now let’s consider a case in which neither \(g\) nor \(f\) generates the data.
Instead, a third distribution \(h\) does.
Let’s study how accumulated likelihood ratios \(L\) behave when \(h\) governs the data.
A key tool here is called Kullback–Leibler divergence we studied in Statistical Divergence Measures.
In our application, we want to measure how much \(f\) or \(g\) diverges from \(h\)
Two Kullback–Leibler divergences pertinent for us are \(K_f\) and \(K_g\) defined as
Let’s compute the Kullback–Leibler discrepancies using the closed form for Beta densities from Statistical Divergence Measures.
def KL_beta(a1, b1, a2, b2):
"""KL(Beta(a1, b1) ‖ Beta(a2, b2)) in closed form."""
return (betaln(a2, b2) - betaln(a1, b1)
+ (a1 - a2) * digamma(a1) + (b1 - b2) * digamma(b1)
+ (a2 - a1 + b2 - b1) * digamma(a1 + b1))
def compute_KL_h(h_params, F_params=(F_a, F_b), G_params=(G_a, G_b)):
"""
Compute K_f = KL(h, f) and K_g = KL(h, g) for Beta densities h, f, g
"""
Kf = KL_beta(*h_params, *F_params)
Kg = KL_beta(*h_params, *G_params)
return Kf, Kg
22.6.2. A helpful formula#
There is a mathematical relationship between likelihood ratios and KL divergence.
When data is generated by distribution \(h\), the expected log likelihood ratio is:
where \(L_t=\prod_{j=1}^{t}\frac{f(w_j)}{g(w_j)}\) is the likelihood ratio process.
Equation (22.1) tells us that:
When \(K_g < K_f\) (i.e., \(g\) is closer to \(h\) than \(f\) is), the expected log likelihood ratio is negative, so \(L\left(w^t\right) \rightarrow 0\).
When \(K_g > K_f\) (i.e., \(f\) is closer to \(h\) than \(g\) is), the expected log likelihood ratio is positive, so \(L\left(w^t\right) \rightarrow + \infty\).
Let’s verify this using simulation.
In the simulation, we generate multiple paths using Beta distributions \(f\), \(g\), and \(h\), and compute the paths of \(\log(L(w^t))\).
First, we write a function to compute the likelihood ratio process
def compute_likelihood_ratios(sequences, f, g):
"""Compute likelihood ratios and cumulative products."""
l_ratios = f(sequences) / g(sequences)
L_cumulative = np.cumprod(l_ratios, axis=1)
return l_ratios, L_cumulative
We consider three cases: (1) \(h\) is closer to \(f\), (2) \(f\) and \(g\) are equidistant from \(h\), and (3) \(h\) is closer to \(g\).
In case (2) we set \(h = \text{Beta}(2.1035, 1.35)\), for which \(K_f\) and \(K_g\) agree to four decimal places.
Note that
In the first figure, \(\log L(w^t)\) diverges to \(\infty\) because \(K_g > K_f\).
In the second figure, \(K_g = K_f\), so \(\log L(w^t)\) is a random walk with IID increments of mean zero; its paths wander away from zero in both directions, with a spread that grows like \(\sqrt{t}\).
In the last figure, \(\log L(w^t)\) diverges to \(-\infty\) because \(K_g < K_f\).
The black dashed line, \(t \left(D_{KL}(h\|g) - D_{KL}(h\|f)\right)\), passes through the center of the paths, as (22.1) predicts.
These observations align with the theory.
In Heterogeneous Beliefs and Financial Markets, we will see an application of these ideas.
22.7. Hypothesis testing and classification#
This section discusses another application of likelihood ratio processes.
We describe how a statistician can combine frequentist probabilities of type I and type II errors in order to
compute an anticipated frequency of selecting a wrong model based on a sample length \(T\)
compute an anticipated error rate in a classification problem
We consider a situation in which nature generates data by mixing known densities \(f\) and \(g\) with known mixing parameter \(\pi_{-1} \in (0,1)\) so that the random variable \(w\) is drawn from the density
We assume that the statistician knows the densities \(f\) and \(g\) and also the mixing parameter \(\pi_{-1}\).
Below, we’ll set \(\pi_{-1} = .5\), although much of the analysis would follow through with other settings of \(\pi_{-1} \in (0,1)\).
We assume that \(f\) and \(g\) both put positive probabilities on the same intervals of possible realizations of the random variable \(W\).
In the simulations below, we specify that \(f\) is a \(\text{Beta}(1, 1)\) distribution and that \(g\) is \(\text{Beta}(3, 1.2)\) distribution.
We consider two alternative timing protocols.
Timing protocol 1 is for the model selection problem
Timing protocol 2 is for the individual classification problem
Timing Protocol 1: Nature flips a coin only once at time \(t=-1\) and with probability \(\pi_{-1}\) generates a sequence \(\{w_t\}_{t=1}^T\) of IID draws from \(f\) and with probability \(1-\pi_{-1}\) generates a sequence \(\{w_t\}_{t=1}^T\) of IID draws from \(g\).
Timing Protocol 2: Nature flips a coin often. At each time \(t \geq 1\), nature flips a coin and with probability \(\pi_{-1}\) draws \(w_t\) from \(f\) and with probability \(1-\pi_{-1}\) draws \(w_t\) from \(g\).
Here is Python code that we’ll use to implement timing protocol 1 and 2
def protocol_1(π_minus_1, T, rng, N=1000, F_params=(1, 1), G_params=(3, 1.2)):
"""
Simulate Protocol 1: Nature decides once at t=-1 which model to use.
"""
F_a, F_b = F_params
G_a, G_b = G_params
# Single coin flip for the true model
true_models_F = rng.random(N) < π_minus_1
sequences = np.empty((N, T))
n_f = np.sum(true_models_F)
n_g = N - n_f
if n_f > 0:
sequences[true_models_F, :] = rng.beta(F_a, F_b, (n_f, T))
if n_g > 0:
sequences[~true_models_F, :] = rng.beta(G_a, G_b, (n_g, T))
return sequences, true_models_F
def protocol_2(π_minus_1, T, rng, N=1000, F_params=(1, 1), G_params=(3, 1.2)):
"""
Simulate Protocol 2: Nature decides at each time step which model to use.
"""
F_a, F_b = F_params
G_a, G_b = G_params
# Coin flips for each time step
true_models_F = rng.random((N, T)) < π_minus_1
sequences = np.empty((N, T))
n_f = np.sum(true_models_F)
n_g = N * T - n_f
if n_f > 0:
sequences[true_models_F] = rng.beta(F_a, F_b, n_f)
if n_g > 0:
sequences[~true_models_F] = rng.beta(G_a, G_b, n_g)
return sequences, true_models_F
Remark: Under timing protocol 2, the \(\{w_t\}_{t=1}^T\) is a sequence of IID draws from \(h(w)\). Under timing protocol 1, the \(\{w_t\}_{t=1}^T\) is not IID. It is conditionally IID – meaning that with probability \(\pi_{-1}\) it is a sequence of IID draws from \(f(w)\) and with probability \(1-\pi_{-1}\) it is a sequence of IID draws from \(g(w)\). For more about this, see this lecture about exchangeability.
Likelihood Ratio Processes and Bayesian Learning studies a Bayesian who knows that timing protocol 1 governs the data and learns which density nature chose.
Incorrect Models studies a Bayesian who wrongly believes timing protocol 1 when timing protocol 2 generates the data.
We again deploy a likelihood ratio process with time \(t\) component being the likelihood ratio
The likelihood ratio process for sequence \(\left\{ w_{t}\right\} _{t=1}^{\infty}\) is
For shorthand we’ll write \(L_t = L(w^t)\).
22.7.1. Model selection mistake probability#
We first study a problem that assumes timing protocol 1.
Consider a decision maker who wants to know whether model \(f\) or model \(g\) governs a data set of length \(T\) observations.
The decision maker has observed a sequence \(\{w_t\}_{t=1}^T\).
On the basis of that observed sequence, a likelihood ratio test selects model \(f\) when \(L_T \geq 1 \) and model \(g\) when \(L_T < 1\).
When model \(f\) generates the data, the probability that the likelihood ratio test selects the wrong model is
When model \(g\) generates the data, the probability that the likelihood ratio test selects the wrong model is
We can construct a probability that the likelihood ratio selects the wrong model by assigning a Bayesian prior probability of \(\pi_{-1} = .5\) that nature selects model \(f\) then averaging \(p_f\) and \(p_g\) to form the prior (ex ante) probability of a detection error equal to
Now let’s simulate timing protocol 1 and compute the error probabilities
def compute_protocol_1_errors(π_minus_1, T_max, N_simulations, f_func, g_func,
rng, F_params=(1, 1), G_params=(3, 1.2)):
"""
Compute error probabilities for Protocol 1.
"""
sequences, true_models = protocol_1(
π_minus_1, T_max, rng, N_simulations, F_params, G_params)
l_ratios, L_cumulative = compute_likelihood_ratios(sequences,
f_func, g_func)
T_range = np.arange(1, T_max + 1)
mask_f = true_models
mask_g = ~true_models
L_f = L_cumulative[mask_f, :]
L_g = L_cumulative[mask_g, :]
α_T = np.mean(L_f < 1, axis=0)
β_T = np.mean(L_g >= 1, axis=0)
error_prob = 0.5 * (α_T + β_T)
return {
'T_range': T_range,
'alpha': α_T,
'beta': β_T,
'error_prob': error_prob,
'L_cumulative': L_cumulative,
'true_models': true_models
}
The following code visualizes the error probabilities for timing protocol 1
At T=30:
α_30 = 0.0046
β_30 = 0.0036
Model selection error probability = 0.0041
Notice how the model selection error probability approaches zero as \(T\) grows.
22.7.2. Classification#
We now consider a problem that assumes timing protocol 2.
A decision maker wants to classify components of an observed sequence \(\{w_t\}_{t=1}^T\) as having been drawn from either \(f\) or \(g\).
The decision maker uses the following classification rule:
Under this rule, the expected misclassification rate is
where \(\tilde \alpha_t = {\rm Prob}(l_t < 1 \mid f)\) and \(\tilde \beta_t = {\rm Prob}(l_t \geq 1 \mid g)\).
Now let’s write some code to simulate it
def compute_protocol_2_errors(π_minus_1, T_max, N_simulations, f_func, g_func,
rng, F_params=(1, 1), G_params=(3, 1.2)):
"""
Compute error probabilities for Protocol 2.
"""
sequences, true_models = protocol_2(π_minus_1,
T_max, rng, N_simulations, F_params, G_params)
l_ratios, _ = compute_likelihood_ratios(sequences, f_func, g_func)
T_range = np.arange(1, T_max + 1)
accuracy = np.empty(T_max)
for t in range(T_max):
predictions = (l_ratios[:, t] >= 1)
actual = true_models[:, t]
accuracy[t] = np.mean(predictions == actual)
return {
'T_range': T_range,
'accuracy': accuracy,
'l_ratios': l_ratios,
'true_models': true_models
}
Since for each \(t\), the decision boundary is the same, the decision boundary can be computed as
root = brentq(lambda w: f(w) / g(w) - 1, 0.001, 0.99)
root_2 = brentq(lambda w: f(w) / g(w) - 1, 0.99, 1 - 1e-12)
print(f"f/g = 1 at w = {root:.4f} and at w = {root_2:.5f}")
f/g = 1 at w = 0.5241 and at w = 0.99925
The likelihood ratio \(f/g\) equals \(1\) at two points, because \(g\) vanishes at both ends of \((0,1)\).
The rule classifies \(w_t\) as coming from \(g\) when \(w_t\) lies between the two roots and as coming from \(f\) otherwise.
The second root lies so close to \(1\) that the interval to its right has probability below \(0.001\) under both \(f\) and \(g\).
We can plot the distributions of \(f\) and \(g\) and the decision boundary
To the left of the green vertical line \(g < f\), so \(l_t > 1\); therefore a \(w_t\) that falls to the left of the green line is classified as a type \(f\) individual.
The shaded red area equals \(\tilde \beta_t\) – the probability of classifying someone as a type \(f\) individual when it is really a type \(g\) individual.
To the right of the green vertical line and to the left of the second root \(w = 0.99925\), \(g > f\), so \(l_t < 1\); therefore a \(w_t\) that falls in this interval is classified as a type \(g\) individual.
The shaded blue area equals \(\tilde \alpha_t\) – the probability of classifying someone as a type \(g\) individual when it is really a type \(f\) individual.
This gives us clues about how to compute the theoretical classification error probability
# tilde α_t = P(l_t < 1 | f) = P(root < w < root_2 | f)
α_theory = beta_dist.cdf(root_2, F_a, F_b) - beta_dist.cdf(root, F_a, F_b)
# tilde β_t = P(l_t >= 1 | g) = P(w <= root or w >= root_2 | g)
β_theory = beta_dist.cdf(root, G_a, G_b) + beta_dist.sf(root_2, G_a, G_b)
theory_error = 0.5 * (α_theory + β_theory)
print(f"theoretical tilde α_t = {α_theory:.4f}")
print(f"theoretical tilde β_t = {β_theory:.4f}")
print(f"theoretical classification error probability = {theory_error:.4f}")
theoretical tilde α_t = 0.4752
theoretical tilde β_t = 0.1836
theoretical classification error probability = 0.3294
Now we simulate timing protocol 2 and compute the classification error probability.
In the next cell, we also compare the theoretical classification accuracy to the empirical classification accuracy
def analyze_protocol_2(π_minus_1, T_max, N_simulations, f_func, g_func, rng,
theory_error=None, F_params=(1, 1), G_params=(3, 1.2)):
"""Analyze Protocol 2."""
result = compute_protocol_2_errors(π_minus_1, T_max, N_simulations,
f_func, g_func, rng, F_params, G_params)
# Plot results
plt.figure(figsize=(10, 6))
plt.plot(result['T_range'], result['accuracy'],
'b-', linewidth=2, label='empirical accuracy')
if theory_error is not None:
plt.axhline(1 - theory_error, color='r', linestyle='--',
label=f'theoretical accuracy = {1 - theory_error:.4f}')
plt.xlabel('$t$')
plt.ylabel('accuracy')
plt.legend()
plt.ylim(0.5, 1.0)
plt.show()
return result
# Analyze Protocol 2
result_p2 = analyze_protocol_2(π_minus_1, T_max, N_simulations, f, g, rng,
theory_error, (F_a, F_b), (G_a, G_b))
Let’s watch decisions made by the two timing protocols as more and more observations accrue.
def compare_protocols(result1, result2):
"""Compare results from both protocols."""
plt.figure(figsize=(10, 6))
plt.plot(result1['T_range'], result1['error_prob'], linewidth=2,
label='Protocol 1 (Model Selection)')
plt.plot(result2['T_range'], 1 - result2['accuracy'],
linestyle='--', linewidth=2,
label='Protocol 2 (classification)')
plt.xlabel('$T$')
plt.ylabel('error probability')
plt.legend()
plt.show()
compare_protocols(result_p1, result_p2)
From the figure above, we can see:
For both timing protocols, the error probability starts at the same level, subject to a little randomness.
For timing protocol 1, the error probability decreases as the sample size increases because we are making just one decision – i.e., selecting whether \(f\) or \(g\) governs all individuals. More data provides better evidence.
For timing protocol 2, the error probability remains constant because we are making many decisions – one classification decision for each observation.
Remark: Think about how laws of large numbers are applied to compute error probabilities for the model selection problem and the classification problem.
22.7.3. Error probability and divergence measures#
A plausible guess is that the ability of a likelihood ratio to distinguish distributions \(f\) and \(g\) depends on how “different” they are.
We have learnt some measures of “difference” between distributions in Statistical Divergence Measures.
Let’s now study two more measures of “difference” between distributions that are useful in the context of model selection and classification.
Recall that Chernoff entropy between probability densities \(f\) and \(g\) is defined as:
Let \(f_T(w^T) = \prod_{t=1}^T f(w_t)\) and \(g_T(w^T) = \prod_{t=1}^T g(w_t)\) be the joint densities of \(w^T\) under the two models.
The likelihood ratio test selects \(f\) when \(f_T \geq g_T\) and \(g\) when \(f_T < g_T\).
Its model selection error probability (22.2) is therefore
For \(x, y \geq 0\) and \(\phi \in [0,1]\), \(\min\{x, y\} \leq x^\phi y^{1-\phi}\).
Because the \(w_t\)’s are IID,
for every \(\phi \in [0,1]\).
Minimizing the right side over \(\phi\) gives an upper bound on model selection error probability:
The weaker bound \(e^{-C(f,g)T}\) also holds.
Chernoff [1952] showed that the bound (22.4) delivers the correct exponential rate:
The approach to the limit is slow, because the error probability contains a factor that decays polynomially in \(T\), which the bound omits.
For Beta densities, the Chernoff coefficient \(\rho(\phi)\) has a closed form that we derived in Statistical Divergence Measures.
Let’s compute Chernoff entropy with some Python code
def log_chernoff_coef_beta(ϕ, a1, b1, a2, b2):
"""log ρ(ϕ) for Beta(a1, b1) and Beta(a2, b2) in closed form."""
return (betaln(ϕ * a1 + (1 - ϕ) * a2, ϕ * b1 + (1 - ϕ) * b2)
- ϕ * betaln(a1, b1) - (1 - ϕ) * betaln(a2, b2))
def chernoff_beta(a1, b1, a2, b2):
"""Chernoff entropy C(f, g) and minimizing ϕ for two Beta densities."""
res = minimize_scalar(log_chernoff_coef_beta, bounds=(0, 1),
args=(a1, b1, a2, b2), method='bounded')
return -res.fun, res.x
C_fg, ϕ_optimal = chernoff_beta(F_a, F_b, G_a, G_b)
print(f"Chernoff entropy C(f,g) = {C_fg:.4f}")
print(f"Optimal ϕ = {ϕ_optimal:.4f}")
Chernoff entropy C(f,g) = 0.1212
Optimal ϕ = 0.5969
Now let’s compare the model selection error probability that we simulated under timing protocol 1 with the bounds \(\frac{1}{2} e^{-C(f,g)T}\) and \(e^{-C(f,g)T}\).
The right panel plots \(-\frac{1}{T} \log\) of the simulated error probability against \(T\).
T_range = np.arange(1, T_max+1)
chernoff_bound = 0.5 * np.exp(-C_fg * T_range)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
ax1.semilogy(T_range, 2 * chernoff_bound, 'r--', linewidth=2,
label=f'$e^{{-C(f,g)T}}$')
ax1.semilogy(T_range, chernoff_bound, 'r-', linewidth=2,
label=r'$\frac{1}{2} e^{-C(f,g)T}$')
ax1.semilogy(T_range, result_p1['error_prob'], 'b-', linewidth=2,
label='Model selection error probability')
ax1.set_xlabel('T')
ax1.set_ylabel('error probability (log scale)')
ax1.legend()
ax2.plot(T_range, -np.log(result_p1['error_prob']) / T_range, 'b-',
linewidth=2, label=r'$-\frac{1}{T}\log \frac{1}{2}(\alpha_T + \beta_T)$')
ax2.axhline(C_fg, color='r', linestyle='-', linewidth=2, label='$C(f,g)$')
ax2.set_xlabel('T')
ax2.set_ylim(0, None)
ax2.legend()
plt.tight_layout()
plt.show()
The simulated error probability lies below \(\frac{1}{2} e^{-C(f,g)T}\), as (22.4) requires.
Its negative log divided by \(T\) declines toward \(C(f,g)\) slowly.
In Statistical Divergence Measures, we also studied Jensen-Shannon divergence as a symmetric measure of distance between distributions.
We can use Jensen-Shannon divergence to measure the distance between distributions \(f\) and \(g\) and compute how it covaries with the model selection error probability.
We also compute Jensen-Shannon divergence numerically with some Python code.
We integrate over all of \((0,1)\), splitting the interval at \(1/2\), because U-shaped Beta densities diverge at the endpoints.
def compute_JS(f, g):
"""
Compute Jensen-Shannon divergence
"""
def integrand(w):
fw, gw = f(w), g(w)
mw = 0.5 * (fw + gw)
return 0.5 * fw * np.log(fw / mw) + 0.5 * gw * np.log(gw / mw)
left, _ = quad(integrand, 0, 0.5, limit=200)
right, _ = quad(integrand, 0.5, 1, limit=200)
return left + right
Now let’s return to our guess that the error probability at large sample sizes is related to the Chernoff entropy between two distributions.
We verify this by computing the correlation between the log of the error probability at \(T=50\) under Timing Protocol 1 and the divergence measures.
We use the following pairs of Beta distributions for \(f\) and \(g\) as test cases
distribution_pairs = [
# (f_params, g_params)
((1, 1), (0.1, 0.2)),
((1, 1), (0.3, 0.3)),
((1, 1), (0.3, 0.4)),
((1, 1), (0.5, 0.5)),
((1, 1), (0.7, 0.6)),
((1, 1), (0.9, 0.8)),
((1, 1), (1.1, 1.05)),
((1, 1), (1.2, 1.1)),
((1, 1), (1.5, 1.2)),
((1, 1), (2, 1.5)),
((1, 1), (2.5, 1.8)),
((1, 1), (3, 1.2)),
((1, 1), (4, 1)),
((1, 1), (5, 1))
]
For the most distinct pair, the Chernoff bound (22.4) puts the error probability at \(T = 50\) below \(\frac{1}{2} e^{-0.443 \times 50} \approx 1.2 \times 10^{-10}\).
A simulation of a few thousand sequences from \(f\) and from \(g\) typically records no errors at all for several of these pairs.
We compute the error probabilities instead by importance sampling, a technique that Mean of a Likelihood Ratio Process describes.
We draw \(w^T\) from the tilted density
where \(\phi^*\) attains the minimum in the definition of \(C(f,g)\).
For Beta densities \(f\) and \(g\), \(q\) is again a Beta density, as Statistical Divergence Measures shows.
Since \(f/q = \rho(\phi^*)\, \ell^{1-\phi^*}\) and \(g/q = \rho(\phi^*)\, \ell^{-\phi^*}\),
where \(q_T(w^T) = \prod_{t=1}^T q(w_t)\).
The random variable inside the expectation lies between \(0\) and \(1\).
That fact delivers the bound (22.4) once more and gives the sample average of \(\min\{L_T^{1-\phi^*}, L_T^{-\phi^*}\}\) a small variance.
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 error_prob_tilted(f_params, g_params, T, N, rng):
"""
Model selection error probability (α_T + β_T) / 2 by
importance sampling from the Chernoff-tilted density.
"""
(f_a, f_b), (g_a, g_b) = f_params, g_params
C, ϕ = chernoff_beta(f_a, f_b, g_a, g_b)
q_a, q_b = ϕ * f_a + (1 - ϕ) * g_a, ϕ * f_b + (1 - ϕ) * g_b
# Keep draws off the boundary, where the log densities are infinite
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)
weights = np.exp(np.minimum((1 - ϕ) * log_L, -ϕ * log_L))
return 0.5 * np.exp(-C * T) * np.mean(weights)
As a check, we apply the importance sampling estimator to our \(f\) and \(g\) at \(T = 30\) and compare it with the simulation of timing protocol 1 above.
is_estimate = error_prob_tilted((F_a, F_b), (G_a, G_b), T_max, 5000, rng)
print(f"Importance sampling estimate at T={T_max}: {is_estimate:.5f}")
print(f"Simulated protocol 1 error at T={T_max}: "
f"{result_p1['error_prob'][-1]:.5f}")
Importance sampling estimate at T=30: 0.00355
Simulated protocol 1 error at T=30: 0.00410
Now we compute the divergence measures and the error probability at \(T = 50\) for every pair
# Parameters for simulation
T_large = 50
N_sims = 5000
# Initialize arrays
n_pairs = len(distribution_pairs)
kl_fg_vals = np.zeros(n_pairs)
kl_gf_vals = np.zeros(n_pairs)
js_vals = np.zeros(n_pairs)
chernoff_vals = np.zeros(n_pairs)
error_probs = np.zeros(n_pairs)
pair_names = []
for i, ((f_a, f_b), (g_a, g_b)) in enumerate(distribution_pairs):
# Create density functions
f_pair = jit(lambda x, a=f_a, b=f_b: p(x, a, b))
g_pair = jit(lambda x, a=g_a, b=g_b: p(x, a, b))
# Compute divergence measures
kl_fg_vals[i] = KL_beta(f_a, f_b, g_a, g_b)
kl_gf_vals[i] = KL_beta(g_a, g_b, f_a, f_b)
js_vals[i] = compute_JS(f_pair, g_pair)
chernoff_vals[i], _ = chernoff_beta(f_a, f_b, g_a, g_b)
# Error probability by importance sampling
error_probs[i] = error_prob_tilted((f_a, f_b), (g_a, g_b),
T_large, N_sims, rng)
pair_names.append(f"Beta({f_a},{f_b}) and Beta({g_a},{g_b})")
cor_data = {
'kl_fg': kl_fg_vals,
'kl_gf': kl_gf_vals,
'js': js_vals,
'chernoff': chernoff_vals,
'error_prob': error_probs,
'names': pair_names,
'T': T_large}
Now let’s visualize the correlations.
The left panel plots the log error probability against \(-C(f,g)T\), together with the Chernoff bound \(\log \frac{1}{2} - C(f,g)T\).
All fourteen points lie below the Chernoff bound.
Evidently, Chernoff entropy and Jensen-Shannon divergence each covary tightly with the model selection error probability.
The vertical distance between the points and the Chernoff bound is the log of \(E_q\left[\min\left\{L_T^{1-\phi^*}, L_T^{-\phi^*}\right\}\right]\), which varies across pairs much less than \(C(f,g) T\) does.
We’ll encounter related ideas in A Problem that Stumped Milton Friedman very soon.
22.8. Markov chains#
Let’s now look at a likelihood ratio process for a sequence of random variables that is not independently and identically distributed.
Here we assume that the sequence is generated by a Markov chain on a finite state space.
We consider two \(n\)-state irreducible and aperiodic Markov chain models on the same state space \(\{1, 2, \ldots, n\}\) with positive transition matrices \(P^{(f)}\), \(P^{(g)}\) and initial distributions \(\pi_0^{(f)}\), \(\pi_0^{(g)}\).
We assume that nature samples from chain \(f\).
For a sample path \((x_0, x_1, \ldots, x_T)\), let \(N_{ij}\) count transitions from state \(i\) to \(j\).
The likelihood process under model \(m \in \{f, g\}\) is
Hence,
The log-likelihood ratio is
22.8.1. KL divergence rate#
By the ergodic theorem for irreducible Markov chains, we have
where \(\boldsymbol{\pi}^{(f)}\) is the stationary distribution satisfying \(\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}\).
Therefore,
Taking the limit as \(T \to \infty\), we have:
The first term: \(\frac{1}{T}\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} \to 0\)
The second term: \(\frac{1}{T}\sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} \xrightarrow{a.s.} \sum_{i,j}\pi_i^{(f)}P_{ij}^{(f)}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}}\)
Define the KL divergence rate as
where \(KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})\) is the row-wise KL divergence.
By the ergodic theorem, we have
Taking expectations and using the dominated convergence theorem, we obtain
Here we invite readers to pause and compare this result with (22.1).
Likelihood Processes For VAR Models extends these calculations to vector autoregressions.
Let’s confirm this in the simulation below.
22.8.2. Simulations#
Let’s implement simulations to illustrate these concepts with a three-state Markov chain.
We start by writing functions to compute the stationary distribution and the KL divergence rate for Markov chain models.
Now let’s create an example with two different 3-state Markov chains.
We are now ready to simulate paths and visualize how likelihood ratios evolve.
We draw \(x_0\) from the stationary distribution \(\pi^{(f)}\) and set \(\pi_0^{(f)} = \pi^{(f)}\) and \(\pi_0^{(g)} = \pi^{(g)}\).
The chain is then stationary under \(f\), so \(E_f N_{ij} = T \pi_i^{(f)} P_{ij}^{(f)}\) and (22.5) implies
for every \(T\).
We verify this formula by plotting both the empirical average and the line predicted by the theory
# Define example Markov chain transition matrices
P_f = np.array([[0.7, 0.2, 0.1],
[0.3, 0.5, 0.2],
[0.1, 0.3, 0.6]])
P_g = np.array([[0.5, 0.3, 0.2],
[0.2, 0.6, 0.2],
[0.2, 0.2, 0.6]])
markov_results = analyze_markov_chains(P_f, P_g, rng)
Stationary distribution (f): [0.41176471 0.32352941 0.26470588]
Stationary distribution (g): [0.28571429 0.38095238 0.33333333]
KL divergence rate h(f, g): 0.0588
KL divergence rate h(g, f): 0.0563
KL(π_f, π_g): 0.0366
22.10. Exercises#
Exercise 22.1
Consider the setting where nature generates data from a third density \(h\).
Let \(\{w_t\}_{t=1}^T\) be IID draws from \(h\), and let \(L_t = L(w^t)\) be the likelihood ratio process defined as in the lecture.
Assume that \(K_f\) and \(K_g\) are finite and that \(E_h |\log f(W)| < \infty\) and \(E_h |\log g(W)| < \infty\).
Show that:
Hint: Start by expressing \(\log L_t\) as a sum of \(\log \ell(w_i)\) terms and compare with the definition of \(K_f\) and \(K_g\).
Solution
Since \(w_1, \ldots, w_t\) are IID draws from \(h\), we can write
Taking the expectation under \(h\)
Since the \(w_i\) are identically distributed
where \(w \sim h\).
Therefore
Now, from the definition of Kullback-Leibler divergence
This gives us
Similarly
Substituting back
Exercise 22.2
Building on Exercise 22.1, use the result to explain what happens to \(L_t\) as \(t \to \infty\) in the following cases:
When \(K_g > K_f\) (i.e., \(f\) is “closer” to \(h\) than \(g\) is)
When \(K_g < K_f\) (i.e., \(g\) is “closer” to \(h\) than \(f\) is)
Relate your answer to the simulation results shown in this section.
Solution
From Exercise 22.1, we know that:
Case 1: When \(K_g > K_f\)
Here, \(f\) is “closer” to \(h\) than \(g\) is. Since \(K_g - K_f > 0\)
By the strong law of large numbers, \(\frac{1}{t} \log L_t \to K_g - K_f > 0\) almost surely.
Therefore \(L_t \to +\infty\) almost surely.
Case 2: When \(K_g < K_f\)
Here, \(g\) is “closer” to \(h\) than \(f\) is. Since \(K_g - K_f < 0\)
Therefore by similar reasoning \(L_t \to 0\) almost surely.
Exercise 22.3
This exercise extends the Chernoff bound (22.4) to the Markov chains of this section.
Let \(x_0\) be drawn from \(\pi_0^{(f)}\) under model \(f\) and from \(\pi_0^{(g)}\) under model \(g\).
The likelihood ratio test selects \(f\) when \(L_T^{(f)}/L_T^{(g)} \geq 1\) and \(g\) otherwise; let \(\alpha_T\) and \(\beta_T\) be its two error probabilities.
For \(\phi \in [0,1]\), let \(Q_\phi\) be the \(n \times n\) matrix with entries \(Q_\phi(i,j) = \bigl(P^{(f)}_{ij}\bigr)^\phi \bigl(P^{(g)}_{ij}\bigr)^{1-\phi}\), and let \(\lambda(\phi)\) be its largest eigenvalue (the Perron root).
Show that
\[ \frac{1}{2}(\alpha_T + \beta_T) \leq \frac{1}{2} v_\phi^\top Q_\phi^T \mathbf{1} , \]where \(v_\phi(i) = \bigl(\pi^{(f)}_{0,i}\bigr)^\phi \bigl(\pi^{(g)}_{0,i}\bigr)^{1-\phi}\) and \(\mathbf{1}\) is a vector of ones.
Use the Perron–Frobenius theorem to conclude that the error probability decays at least at the exponential rate \(C_M = -\log \min_{\phi \in [0,1]} \lambda(\phi)\).
Compute \(C_M\) and the minimizing \(\phi\) for the matrices \(P^{(f)}\) and \(P^{(g)}\) in the lecture, and compare \(C_M\) with the KL divergence rates \(h_{KL}(f,g)\) and \(h_{KL}(g,f)\).
Set \(\pi_0^{(f)} = \pi^{(f)}\) and \(\pi_0^{(g)} = \pi^{(g)}\). Use
compute_markov_selection_errorto simulate the error probability for several values of \(T\), and plot it together with the bound in part 1.
Solution
Part 1. Let \(p_f(x^T) = \pi^{(f)}_{0,x_0} \prod_{t=1}^T P^{(f)}_{x_{t-1} x_t}\) and define \(p_g(x^T)\) in the same way.
As in the IID case, the error probability equals \(\frac{1}{2} \sum_{x^T} \min\{p_f(x^T), p_g(x^T)\}\).
Since \(\min\{x, y\} \leq x^\phi y^{1-\phi}\),
Part 2. The matrix \(Q_\phi\) has positive entries.
By the Perron–Frobenius theorem, \(Q_\phi^T = \lambda(\phi)^T \bigl(r_\phi \ell_\phi^\top + o(1)\bigr)\) as \(T \rightarrow \infty\), where \(r_\phi\) and \(\ell_\phi\) are positive right and left eigenvectors normalized so that \(\ell_\phi^\top r_\phi = 1\).
Hence \(\frac{1}{T} \log \bigl(v_\phi^\top Q_\phi^T \mathbf{1}\bigr) \rightarrow \log \lambda(\phi)\), and the error probability decays at least at rate \(-\log \lambda(\phi)\) for every \(\phi\).
The best rate comes from the \(\phi\) that minimizes \(\lambda(\phi)\).
Because \(P^{(f)}\) and \(P^{(g)}\) are stochastic matrices, \(\lambda(0) = \lambda(1) = 1\).
Large deviation theory for Markov chains implies that \(C_M\) is the exact rate [Dembo and Zeitouni, 1998].
Part 3.
def perron_root(ϕ, P_f, P_g):
"""Largest eigenvalue of the matrix with entries P_f^ϕ P_g^(1-ϕ)."""
Q = P_f**ϕ * P_g**(1 - ϕ)
return np.max(np.abs(np.linalg.eigvals(Q)))
res = minimize_scalar(perron_root, bounds=(0, 1), args=(P_f, P_g),
method='bounded')
C_M, ϕ_M = -np.log(res.fun), res.x
print(f"C_M = {C_M:.5f} at ϕ = {ϕ_M:.4f}")
print(f"h_KL(f, g) = {markov_results['kl_rate_fg']:.4f}, "
f"h_KL(g, f) = {markov_results['kl_rate_gf']:.4f}")
C_M = 0.01427 at ϕ = 0.5055
h_KL(f, g) = 0.0588, h_KL(g, f) = 0.0563
The Chernoff rate \(C_M = 0.0143\) is about a quarter of either KL divergence rate.
The error probability therefore falls much more slowly than \(\log L_T^{(f)}/L_T^{(g)}\) drifts.
Part 4.
π_f = markov_results['stationary_f']
π_g = markov_results['stationary_g']
T_values = np.array([25, 50, 100, 150, 200, 250, 300])
errors = compute_markov_selection_error(T_values, P_f, P_g, π_f, π_g,
rng, N_sim=4000)
Q_M = P_f**ϕ_M * P_g**(1 - ϕ_M)
v_M = π_f**ϕ_M * π_g**(1 - ϕ_M)
bound = [0.5 * v_M @ np.linalg.matrix_power(Q_M, T) @ np.ones(3)
for T in T_values]
plt.figure(figsize=(10, 6))
plt.semilogy(T_values, errors, 'bo-', label='simulated error probability')
plt.semilogy(T_values, bound, 'r-',
label=r'bound $\frac{1}{2} v_\phi^\top Q_\phi^T \mathbf{1}$')
plt.semilogy(T_values, 0.5 * np.exp(-C_M * T_values), 'r--',
label=r'$\frac{1}{2} e^{-C_M T}$')
plt.xlabel('T')
plt.ylabel('error probability (log scale)')
plt.legend()
plt.show()
The bound from part 1 nearly coincides with \(\frac{1}{2} e^{-C_M T}\).
The simulated error probabilities lie below the bound.
On the log scale they decline roughly in parallel with the bound, a little faster because of a polynomial factor that the bound omits.
Exercise 22.4
Let \(f = \text{Beta}(1,1)\) and \(g = \text{Beta}(G_a, G_b)\).
Show that
\[ E\left[\ell(w_t) \mid q = f\right] = \int_0^1 \frac{f(w)^2}{g(w)} dw = B(G_a, G_b) \, B(2 - G_a, 2 - G_b) \]when \(G_a < 2\) and \(G_b < 2\), and that \(E\left[\ell(w_t) \mid q = f\right] = +\infty\) otherwise.
Compute \(E\left[\ell(w_t) \mid q = f\right]\) for \(g = \text{Beta}(1.5, 1.2)\) and for the lecture’s \(g = \text{Beta}(3, 1.2)\).
For each \(g\), draw \(w_1, w_2, \ldots, w_N\) from \(f\) and plot the running sample average \(\frac{1}{n} \sum_{i=1}^n \ell(w_i)\) against \(n\). Explain the difference between the two plots.
Solution
Part 1. Since \(f = 1\) on \((0,1)\),
The integral converges if and only if \(1 - G_a > -1\) and \(1 - G_b > -1\), that is, if and only if \(G_a < 2\) and \(G_b < 2\).
When it converges, it equals \(B(2 - G_a, 2 - G_b)\).
Part 2.
def mean_ell_under_f(G_a, G_b):
"""E[ℓ(w) | q=f] for f = Beta(1, 1) and g = Beta(G_a, G_b)."""
if G_a < 2 and G_b < 2:
return np.exp(betaln(G_a, G_b) + betaln(2 - G_a, 2 - G_b))
return np.inf
for G_params in [(1.5, 1.2), (3, 1.2)]:
print(f"g = Beta{G_params}: E[ℓ | f] = {mean_ell_under_f(*G_params):.4f}")
g = Beta(1.5, 1.2): E[ℓ | f] = 1.2112
g = Beta(3, 1.2): E[ℓ | f] = inf
Part 3.
N = 1_000_000
n_grid = np.arange(1, N+1)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for ax, (g_a, g_b) in zip(axes, [(1.5, 1.2), (3, 1.2)]):
w = rng.uniform(size=N)
ℓ = np.exp(-log_beta_pdf(w, g_a, g_b))
ax.plot(n_grid, np.cumsum(ℓ) / n_grid, lw=1.5,
label='running sample average')
mean = mean_ell_under_f(g_a, g_b)
if np.isfinite(mean):
ax.axhline(mean, color='k', ls='--', label=r'$E[\ell \mid q=f]$')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('n')
ax.set_title(f'g = Beta({g_a}, {g_b})')
ax.legend()
plt.tight_layout()
plt.show()
For \(g = \text{Beta}(1.5, 1.2)\), the mean of \(\ell\) under \(f\) is finite and the strong law of large numbers drives the running average to \(1.2112\).
For \(g = \text{Beta}(3, 1.2)\), the mean is infinite.
Rare draws near \(w = 0\) produce enormous values of \(\ell = 1/g\), and the running average keeps jumping upward without settling.
For this reason we plotted \(\log L_t\) in the section in which nature draws from f.