From 4d4c572ac7d0da948a2025a9288de702ac2f6ad3 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Sun, 3 Aug 2025 21:53:20 +1000 Subject: [PATCH 1/6] update first draft --- lectures/likelihood_ratio_process.md | 270 ++++++++++++++++++++++++++- 1 file changed, 266 insertions(+), 4 deletions(-) diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index a2037eae6..d1bf32932 100644 --- a/lectures/likelihood_ratio_process.md +++ b/lectures/likelihood_ratio_process.md @@ -4,7 +4,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.2 + jupytext_version: 1.17.1 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -57,6 +57,7 @@ from scipy.optimize import brentq, minimize_scalar from scipy.stats import beta as beta_dist import pandas as pd from IPython.display import display, Math +import quantecon as qe ``` ## Likelihood Ratio Process @@ -1764,6 +1765,270 @@ From the figure above, we can see: **Remark:** Think about how laws of large numbers are applied to compute error probabilities for the model selection problem and the classification problem. +## Special case: Markov chain models + +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)}$. + +In this section, we assume nature chooses $f$. + +For a sample path $(x_0, x_1, \ldots, x_T)$, let $N$ be the transition count matrix where $N_{ij} = n_{ij}$ counts transitions from state $i$ to $j$. + +The likelihood under model $m \in \{f, g\}$ is: + +$$ +L_T^{(m)} = \pi_{0,x_0}^{(m)} \prod_{i=1}^n \prod_{j=1}^n \left(P_{ij}^{(m)}\right)^{n_{ij}} +$$ + +Hence, $\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}n_{ij}\log P_{ij}^{(m)}$. + +In matrix form, we have + +$$ +\log L_T^{(m)} = \log \pi_{0,x_0}^{(m)} + \text{tr}(N^T \log P^{(m)}) +$$ + +where $\log P^{(m)}$ denotes element-wise logarithm and $\text{tr}(\cdot)$ is the trace operator. + +Subtracting two log-likelihoods gives the log-likelihood ratio + +$$ +\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \text{tr}\left(N^T \log \frac{P^{(f)}}{P^{(g)}}\right) +$$ + +where the division is element-wise. + +For an irreducible, aperiodic finite chain, the ergodic theorem ensures that as $T\to\infty$ + +$$ +\frac{n_{ij}}{T} \xrightarrow{a.s.} \pi_i^{(f)}P_{ij}^{(f)}, +$$ + +where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}$. + +### KL divergence rate + +From the log-likelihood ratio formula, taking expectations under model $f$: + +$$ +E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = E_f\left[\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}}\right] + \sum_{i,j} E_f[n_{ij}] \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +$$ + +By the ergodic theorem, $E_f[n_{ij}] = T\pi_i^{(f)}P_{ij}^{(f)} + o(T)$, hence: + +$$ +E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T\sum_{i,j} \pi_i^{(f)}P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} + O(1) +$$ + +Define the **row-wise KL divergence** between the $i$-th rows of the transition matrices + +$$ +D_{KL}(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) := \sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +$$ + +The leading term is the **KL divergence rate** + +$$ +h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} D_{KL}(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) +$$ + +Therefore: + +$$ +E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T \cdot h_{KL}(f, g) + E_f\left[\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}}\right] +$$ + +The initial distribution term is $O(1)$. Thus: + +$$ +\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] \to h_{KL}(f, g) \quad \text{as } T \to \infty +$$ + +We'll confirm this in the simulation below. + +### Simulations + +Let's implement simulations to illustrate these concepts with a three-state Markov chain. + +We start the simulation from the stationary distribution + +```{code-cell} ipython3 +def compute_stationary_dist(P): + """ + Compute stationary distribution of transition matrix P + """ + + eigenvalues, eigenvectors = np.linalg.eig(P.T) + idx = np.argmax(np.abs(eigenvalues)) + stationary = np.real(eigenvectors[:, idx]) + return stationary / stationary.sum() + +def markov_kl_divergence(P_f, P_g, pi_f): + """ + Compute KL divergence rate between two Markov chains + """ + kl_rate = 0.0 + for i in range(len(P_f)): + for j in range(len(P_f)): + if P_f[i, j] > 0 and P_g[i, j] > 0: + kl_rate += pi_f[i] * P_f[i, j] * np.log(P_f[i, j] / P_g[i, j]) + else: + return np.inf + + return kl_rate + +def simulate_markov_chain(P, pi_0, T, N_paths=1000): + """ + Simulate N_paths sample paths from a Markov chain + """ + mc = qe.MarkovChain(P, state_values=None) + + initial_states = np.random.choice(len(P), size=N_paths, p=pi_0) + + paths = np.zeros((N_paths, T+1), dtype=int) + + for i in range(N_paths): + path = mc.simulate(T+1, init=initial_states[i]) + paths[i, :] = path + + return paths + + +def compute_likelihood_ratio_markov(paths, P_f, P_g, π_0_f, π_0_g): + """ + Compute likelihood ratio process for Markov chain paths + """ + N_paths, T_plus_1 = paths.shape + T = T_plus_1 - 1 + L_ratios = np.ones((N_paths, T+1)) + + L_ratios[:, 0] = π_0_f[paths[:, 0]] / π_0_g[paths[:, 0]] + + for t in range(1, T+1): + prev_states = paths[:, t-1] + curr_states = paths[:, t] + + transition_ratios = P_f[prev_states, curr_states] / P_g[prev_states, curr_states] + L_ratios[:, t] = L_ratios[:, t-1] * transition_ratios + + return L_ratios +``` + +Now let's create an example with two different 3-state Markov chains: + +```{code-cell} ipython3 +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]]) + +# Initial distributions +π_0_f = np.array([0.33, 0.33, 0.34]) +π_0_g = np.array([0.33, 0.33, 0.34]) + +# Compute stationary distributions +pi_f = compute_stationary_dist(P_f) +pi_g = compute_stationary_dist(P_g) + +print(f"Stationary distribution (f): {pi_f}") +print(f"Stationary distribution (g): {pi_g}") + +# Compute KL divergence rate +kl_rate_fg = markov_kl_divergence(P_f, P_g, pi_f) +kl_rate_gf = markov_kl_divergence(P_g, P_f, pi_g) + +print(f"\nKL divergence rate h(f, g): {kl_rate_fg:.4f}") +print(f"KL divergence rate h(g, f): {kl_rate_gf:.4f}") +``` + +Let's simulate paths and visualize how likelihood ratios evolve. + +We'll verify $\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = h_{KL}(f, g)$ starting from the stationary distribution by plotting both the empirical average and the theoretical expectation + +```{code-cell} ipython3 +# Simulate paths from model f +T = 100 +N_paths = 1000 +paths_from_f = simulate_markov_chain(P_f, π_0_f, T, N_paths) + +L_ratios_f = compute_likelihood_ratio_markov(paths_from_f, + P_f, P_g, + π_0_f, π_0_g) + +plt.figure(figsize=(10, 6)) + +# Plot individual paths +n_show = 50 +for i in range(n_show): + plt.plot(np.log(L_ratios_f[i, :]), + alpha=0.3, color='blue', lw=0.8) + +# Compute theoretical expectation +theory_line = kl_rate_fg * np.arange(T+1) +plt.plot(theory_line, 'k--', linewidth=2.5, + label=r'$T \times h_{KL}(f,g)$') + +# Compute empirical mean +avg_log_L = np.mean(np.log(L_ratios_f), axis=0) +plt.plot(avg_log_L, 'r-', linewidth=2.5, + label='empirical average', alpha=0.8) + +plt.axhline(y=0, color='gray', + linestyle='--', alpha=0.5) +plt.xlabel(r'$T$') +plt.ylabel(r'$\log L_T$') +plt.title('nature $= f$') +plt.legend() +plt.show() +``` + +Let's examine how the model selection error probability depends on sample size using the same simulation strategy in the previous section + +```{code-cell} ipython3 +def compute_selection_error(T_values, P_f, P_g, π_0_f, π_0_g, N_sim=1000): + """ + Compute model selection error probability for different sample sizes + """ + errors = [] + + for T in T_values: + # Simulate from both models + paths_f = simulate_markov_chain(P_f, π_0_f, T, N_sim//2) + paths_g = simulate_markov_chain(P_g, π_0_g, T, N_sim//2) + + # Compute likelihood ratios + L_f = compute_likelihood_ratio_markov(paths_f, + P_f, P_g, + π_0_f, π_0_g) + L_g = compute_likelihood_ratio_markov(paths_g, + P_f, P_g, + π_0_f, π_0_g) + + # Decision rule: choose f if L_T >= 1 + error_f = np.mean(L_f[:, -1] < 1) # Type I error + error_g = np.mean(L_g[:, -1] >= 1) # Type II error + + total_error = 0.5 * (error_f + error_g) + errors.append(total_error) + + return np.array(errors) + +# Compute error probabilities +T_values = np.arange(10, 201, 10) +errors = compute_selection_error(T_values, + P_f, P_g, + π_0_f, π_0_g) + +# Plot results +plt.figure(figsize=(10, 6)) +plt.plot(T_values, errors, linewidth=2) +plt.xlabel('$T$') +plt.ylabel('model selection error probability') +plt.show() +``` + ## Measuring discrepancies between distributions A plausible guess is that the ability of a likelihood ratio to distinguish distributions $f$ and $g$ depends on how "different" they are. @@ -2405,6 +2670,3 @@ $$ ```{solution-end} ``` - - - From 71a1501414731bee5d207e5a22547721c3aacfb5 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 4 Aug 2025 19:46:01 +1000 Subject: [PATCH 2/6] updates --- lectures/likelihood_ratio_process.md | 96 +++++++++++++--------------- 1 file changed, 46 insertions(+), 50 deletions(-) diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index d1bf32932..e40edbce8 100644 --- a/lectures/likelihood_ratio_process.md +++ b/lectures/likelihood_ratio_process.md @@ -1779,24 +1779,19 @@ $$ L_T^{(m)} = \pi_{0,x_0}^{(m)} \prod_{i=1}^n \prod_{j=1}^n \left(P_{ij}^{(m)}\right)^{n_{ij}} $$ -Hence, $\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}n_{ij}\log P_{ij}^{(m)}$. - -In matrix form, we have +Hence, $$ -\log L_T^{(m)} = \log \pi_{0,x_0}^{(m)} + \text{tr}(N^T \log P^{(m)}) +\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}n_{ij}\log P_{ij}^{(m)} $$ -where $\log P^{(m)}$ denotes element-wise logarithm and $\text{tr}(\cdot)$ is the trace operator. Subtracting two log-likelihoods gives the log-likelihood ratio $$ -\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \text{tr}\left(N^T \log \frac{P^{(f)}}{P^{(g)}}\right) -$$ - -where the division is element-wise. - +\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \sum_{i,j}n_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +$$ (eq:llr_markov) + For an irreducible, aperiodic finite chain, the ergodic theorem ensures that as $T\to\infty$ $$ @@ -1807,49 +1802,51 @@ where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\bolds ### KL divergence rate -From the log-likelihood ratio formula, taking expectations under model $f$: +From {eq}`eq:llr_markov`, taking expectations under model $f$: $$ -E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = E_f\left[\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}}\right] + \sum_{i,j} E_f[n_{ij}] \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = E_f\left[\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}}\right] + \sum_{i,j} E_f[n_{ij}] \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} $$ By the ergodic theorem, $E_f[n_{ij}] = T\pi_i^{(f)}P_{ij}^{(f)} + o(T)$, hence: $$ -E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T\sum_{i,j} \pi_i^{(f)}P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} + O(1) +E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T\sum_{i,j} \pi_i^{(f)}P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} + o(T) $$ Define the **row-wise KL divergence** between the $i$-th rows of the transition matrices $$ -D_{KL}(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) := \sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +KL(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) := \sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} $$ -The leading term is the **KL divergence rate** +Weighed by the stationary distribution, we obtain the **KL divergence rate** $$ -h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} D_{KL}(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) +h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} KL(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) $$ -Therefore: +Therefore, $$ -E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T \cdot h_{KL}(f, g) + E_f\left[\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}}\right] +E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T \cdot h_{KL}(f, g) + o(T) $$ -The initial distribution term is $O(1)$. Thus: +Thus, $$ \frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] \to h_{KL}(f, g) \quad \text{as } T \to \infty $$ -We'll confirm this in the simulation below. +Here we invite readers to pause and compare this result with {eq}`eq:kl_likelihood_link`. + +Let's confirm this in the simulation below. ### Simulations Let's implement simulations to illustrate these concepts with a three-state Markov chain. -We start the simulation from the stationary distribution +We start with writing out functions to compute the stationary distribution and the KL divergence rate for Markov chain models ```{code-cell} ipython3 def compute_stationary_dist(P): @@ -1866,13 +1863,16 @@ def markov_kl_divergence(P_f, P_g, pi_f): """ Compute KL divergence rate between two Markov chains """ - kl_rate = 0.0 - for i in range(len(P_f)): - for j in range(len(P_f)): - if P_f[i, j] > 0 and P_g[i, j] > 0: - kl_rate += pi_f[i] * P_f[i, j] * np.log(P_f[i, j] / P_g[i, j]) - else: - return np.inf + if np.any((P_f > 0) & (P_g == 0)): + return np.inf + + valid_mask = (P_f > 0) & (P_g > 0) + + log_ratios = np.zeros_like(P_f) + log_ratios[valid_mask] = np.log(P_f[valid_mask] / P_g[valid_mask]) + + # Weight by stationary probabilities and sum + kl_rate = np.sum(pi_f[:, np.newaxis] * P_f * log_ratios) return kl_rate @@ -1907,13 +1907,14 @@ def compute_likelihood_ratio_markov(paths, P_f, P_g, π_0_f, π_0_g): prev_states = paths[:, t-1] curr_states = paths[:, t] - transition_ratios = P_f[prev_states, curr_states] / P_g[prev_states, curr_states] + transition_ratios = P_f[prev_states, curr_states] \ + / P_g[prev_states, curr_states] L_ratios[:, t] = L_ratios[:, t-1] * transition_ratios return L_ratios ``` -Now let's create an example with two different 3-state Markov chains: +Now let's create an example with two different 3-state Markov chains ```{code-cell} ipython3 P_f = np.array([[0.7, 0.2, 0.1], @@ -1924,38 +1925,33 @@ P_g = np.array([[0.5, 0.3, 0.2], [0.2, 0.6, 0.2], [0.2, 0.2, 0.6]]) -# Initial distributions -π_0_f = np.array([0.33, 0.33, 0.34]) -π_0_g = np.array([0.33, 0.33, 0.34]) - # Compute stationary distributions -pi_f = compute_stationary_dist(P_f) -pi_g = compute_stationary_dist(P_g) +π_f = compute_stationary_dist(P_f) +π_g = compute_stationary_dist(P_g) -print(f"Stationary distribution (f): {pi_f}") -print(f"Stationary distribution (g): {pi_g}") +print(f"Stationary distribution (f): {π_f}") +print(f"Stationary distribution (g): {π_g}") # Compute KL divergence rate -kl_rate_fg = markov_kl_divergence(P_f, P_g, pi_f) -kl_rate_gf = markov_kl_divergence(P_g, P_f, pi_g) +kl_rate_fg = markov_kl_divergence(P_f, P_g, π_f) +kl_rate_gf = markov_kl_divergence(P_g, P_f, π_g) print(f"\nKL divergence rate h(f, g): {kl_rate_fg:.4f}") print(f"KL divergence rate h(g, f): {kl_rate_gf:.4f}") ``` -Let's simulate paths and visualize how likelihood ratios evolve. +We are now ready to simulate paths and visualize how likelihood ratios evolve. -We'll verify $\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = h_{KL}(f, g)$ starting from the stationary distribution by plotting both the empirical average and the theoretical expectation +We'll verify $\frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = h_{KL}(f, g)$ starting from the stationary distribution by plotting both the empirical average and the line predicted by the theory ```{code-cell} ipython3 -# Simulate paths from model f -T = 100 +T = 500 N_paths = 1000 -paths_from_f = simulate_markov_chain(P_f, π_0_f, T, N_paths) +paths_from_f = simulate_markov_chain(P_f, π_f, T, N_paths) L_ratios_f = compute_likelihood_ratio_markov(paths_from_f, P_f, P_g, - π_0_f, π_0_g) + π_f, π_g) plt.figure(figsize=(10, 6)) @@ -1973,13 +1969,13 @@ plt.plot(theory_line, 'k--', linewidth=2.5, # Compute empirical mean avg_log_L = np.mean(np.log(L_ratios_f), axis=0) plt.plot(avg_log_L, 'r-', linewidth=2.5, - label='empirical average', alpha=0.8) + label='empirical average', alpha=0.5) plt.axhline(y=0, color='gray', linestyle='--', alpha=0.5) plt.xlabel(r'$T$') plt.ylabel(r'$\log L_T$') -plt.title('nature $= f$') +plt.title('nature = $f$') plt.legend() plt.show() ``` @@ -2019,13 +2015,13 @@ def compute_selection_error(T_values, P_f, P_g, π_0_f, π_0_g, N_sim=1000): T_values = np.arange(10, 201, 10) errors = compute_selection_error(T_values, P_f, P_g, - π_0_f, π_0_g) + π_f, π_g) # Plot results plt.figure(figsize=(10, 6)) plt.plot(T_values, errors, linewidth=2) plt.xlabel('$T$') -plt.ylabel('model selection error probability') +plt.ylabel('error probability') plt.show() ``` From cd77186008a26e8ca1b160d18508742aaeb9887d Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 4 Aug 2025 22:35:52 +1000 Subject: [PATCH 3/6] updates --- lectures/likelihood_ratio_process.md | 52 ++++++++++++---------------- 1 file changed, 23 insertions(+), 29 deletions(-) diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index e40edbce8..b8c2b2909 100644 --- a/lectures/likelihood_ratio_process.md +++ b/lectures/likelihood_ratio_process.md @@ -1771,68 +1771,62 @@ Consider two $n$-state irreducible and aperiodic Markov chain models on the same In this section, we assume nature chooses $f$. -For a sample path $(x_0, x_1, \ldots, x_T)$, let $N$ be the transition count matrix where $N_{ij} = n_{ij}$ counts transitions from state $i$ to $j$. +For a sample path $(x_0, x_1, \ldots, x_T)$, let $N_{ij}$ count transitions from state $i$ to $j$. -The likelihood under model $m \in \{f, g\}$ is: +The likelihood under model $m \in \{f, g\}$ is $$ -L_T^{(m)} = \pi_{0,x_0}^{(m)} \prod_{i=1}^n \prod_{j=1}^n \left(P_{ij}^{(m)}\right)^{n_{ij}} +L_T^{(m)} = \pi_{0,x_0}^{(m)} \prod_{i=1}^n \prod_{j=1}^n \left(P_{ij}^{(m)}\right)^{N_{ij}} $$ Hence, $$ -\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}n_{ij}\log P_{ij}^{(m)} +\log L_T^{(m)} =\log\pi_{0,x_0}^{(m)} +\sum_{i,j}N_{ij}\log P_{ij}^{(m)} $$ - -Subtracting two log-likelihoods gives the log-likelihood ratio +The log-likelihood ratio is $$ -\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \sum_{i,j}n_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +\log \frac{L_T^{(f)}}{L_T^{(g)}} = \log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} $$ (eq:llr_markov) - -For an irreducible, aperiodic finite chain, the ergodic theorem ensures that as $T\to\infty$ - -$$ -\frac{n_{ij}}{T} \xrightarrow{a.s.} \pi_i^{(f)}P_{ij}^{(f)}, -$$ - -where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}$. ### KL divergence rate -From {eq}`eq:llr_markov`, taking expectations under model $f$: +By the ergodic theorem for irreducible, aperiodic Markov chains, we have $$ -E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = E_f\left[\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}}\right] + \sum_{i,j} E_f[n_{ij}] \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +\frac{N_{ij}}{T} \xrightarrow{a.s.} \pi_i^{(f)}P_{ij}^{(f)} \quad \text{as } T \to \infty $$ -By the ergodic theorem, $E_f[n_{ij}] = T\pi_i^{(f)}P_{ij}^{(f)} + o(T)$, hence: - -$$ -E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T\sum_{i,j} \pi_i^{(f)}P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} + o(T) -$$ +where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}$. -Define the **row-wise KL divergence** between the $i$-th rows of the transition matrices +Therefore, $$ -KL(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) := \sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} +\frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} = \frac{1}{T}\log \frac{\pi_{0,x_0}^{(f)}}{\pi_{0,x_0}^{(g)}} + \frac{1}{T}\sum_{i,j}N_{ij}\log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}} $$ -Weighed by the stationary distribution, we obtain the **KL divergence rate** +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 $$ -h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} KL(P_{i,\cdot}^{(f)}, P_{i,\cdot}^{(g)}) +h_{KL}(f, g) = \sum_{i=1}^n \pi_i^{(f)} \underbrace{\sum_{j=1}^n P_{ij}^{(f)} \log \frac{P_{ij}^{(f)}}{P_{ij}^{(g)}}}_{=: KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})} $$ -Therefore, +where $KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})$ is the row-wise KL divergence. + + +By the strong law of large numbers for Markov chains, we have $$ -E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] = T \cdot h_{KL}(f, g) + o(T) +\frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} \xrightarrow{a.s.} h_{KL}(f, g) \quad \text{as } T \to \infty $$ -Thus, +Taking expectations and using the dominated convergence theorem, we can get $$ \frac{1}{T}E_f\left[\log \frac{L_T^{(f)}}{L_T^{(g)}}\right] \to h_{KL}(f, g) \quad \text{as } T \to \infty From 7009bf389cd7ee403bace9fffe9efa7b1cf763e6 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Mon, 4 Aug 2025 22:37:57 +1000 Subject: [PATCH 4/6] update --- lectures/likelihood_ratio_process.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index b8c2b2909..c9a7357bb 100644 --- a/lectures/likelihood_ratio_process.md +++ b/lectures/likelihood_ratio_process.md @@ -1820,7 +1820,7 @@ $$ where $KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})$ is the row-wise KL divergence. -By the strong law of large numbers for Markov chains, we have +By the ergodic theorem, we have $$ \frac{1}{T}\log \frac{L_T^{(f)}}{L_T^{(g)}} \xrightarrow{a.s.} h_{KL}(f, g) \quad \text{as } T \to \infty From 7fefc4b5d8caebaee518a207a7b3dc96197eadc4 Mon Sep 17 00:00:00 2001 From: thomassargent30 Date: Tue, 5 Aug 2025 12:07:11 -0600 Subject: [PATCH 5/6] Tom's Aug 5 edits of navy_captain lecture --- lectures/navy_captain.md | 108 +++++++++++++++++++-------------------- 1 file changed, 54 insertions(+), 54 deletions(-) diff --git a/lectures/navy_captain.md b/lectures/navy_captain.md index 56431f556..0846428b6 100644 --- a/lectures/navy_captain.md +++ b/lectures/navy_captain.md @@ -39,6 +39,7 @@ from scipy.optimize import minimize This lecture follows up on ideas presented in the following lectures: * {doc}`A Problem that Stumped Milton Friedman ` +* {doc}`A Bayesian Formulation of Friedman and Wald's Problem ` * {doc}`Exchangeability and Bayesian Updating ` * {doc}`Likelihood Ratio Processes ` @@ -60,9 +61,9 @@ this lecture {doc}`Exchangeability and Bayesian Updating ` and in {doc}`Likelihood Ratio Processes `, which describes the link between Bayesian updating and likelihood ratio processes. -The present lecture uses Python to generate simulations that evaluate expected losses under **frequentist** and **Bayesian** decision rules for an instance of the Navy Captain's decision problem. +The present lecture uses Python to generate simulations that evaluate expected losses under the Neyman-Pearson **frequentist** procedure that the Navy captain questioned and the **Bayesian** decision rule described in {doc}`A Bayesian Formulation of Friedman and Wald's Problem `. -The simulations confirm the Navy Captain's hunch that there is a better rule than the one the Navy had ordered him to use. +The simulations confirm the Navy Captain's hunch that there is a better rule than the Neyman-Pearson likelihood ratio test that the Navy had told him to use. ## Setup @@ -94,9 +95,9 @@ The decision maker pays a cost $c$ for drawing another $z$. We mainly borrow parameters from the quantecon lecture -{doc}`A Problem that Stumped Milton Friedman ` except that we increase both $\bar L_{0}$ + {doc}`A Bayesian Formulation of Friedman and Wald's Problem ` except that we increase both $\bar L_{0}$ and $\bar L_{1}$ from $25$ to $100$ to encourage the -frequentist Navy Captain to take more draws before deciding. +Bayesian decision rule to take more draws before deciding. We set the cost $c$ of taking one more draw at $1.25$. @@ -206,11 +207,7 @@ $f_1$ ## Frequentist Decision Rule -The Navy told the Captain to use a frequentist decision rule. - -In particular, it gave him a decision rule that the Navy had designed by using -frequentist statistical theory to minimize an -expected loss function. +The Navy told the Captain to use a Neyman-Pearson likelihood ratio decision rule. That decision rule is characterized by @@ -228,7 +225,10 @@ The decision rule associated with a sample size $t$ is: is greater than $d$ - decide that $f_1$ is the distribution if the likelihood ratio is less than $d$ -To understand how that rule was engineered, let null and alternative +For our purposes here, we want to compute an expected loss from using this rule, where we borrow +loss parameters $\bar L_1$ and $\bar L_2$ from {doc}`A Bayesian Formulation of Friedman and Wald's Problem `. + +Let null and alternative hypotheses be - null: $H_{0}$: $f=f_{0}$, @@ -269,6 +269,8 @@ To solve for $\bar{V}_{fre}\left(t,d\right)$ numerically, we first simulate sequences of $z$ when either $f_0$ or $f_1$ generates data. +Let's plot empirical distributions, i.e., histograms, associated with $f_0$ and $f_1$. + ```{code-cell} python3 N = 10000 T = 100 @@ -325,7 +327,7 @@ plt.title("Receiver Operating Characteristic Curve") plt.show() ``` -Our frequentist minimizes the expected total loss presented in equation +We can minimizes the expected total loss presented in equation {eq}`val1` by choosing $\left(t,d\right)$. Doing that delivers an expected loss @@ -404,10 +406,8 @@ plt.show() t_optimal = np.argmin(V_fre_arr) + 1 ``` -```{code-cell} python3 -msg = f"The above graph indicates that minimizing over t tells the frequentist to draw {t_optimal} observations and then decide." -print(msg) -``` + +The above graph illustrates how minimizing over $t$ tells the frequentist to draw $t_{\rm optimal}$ observations and then decide. Let’s now change the value of $\pi^{*}$ and watch how the decision rule changes. @@ -461,17 +461,17 @@ plt.show() ## Bayesian Decision Rule In {doc}`A Problem that Stumped Milton Friedman `, -we learned how Abraham Wald confirmed the Navy -Captain’s hunch that there is a better decision rule. +we learned how Abraham Wald confirmed the Navy Captain’s hunch that there is a better decision rule. -We presented a Bayesian procedure that instructed the Captain to makes -decisions by comparing his current Bayesian posterior probability -$\pi$ with two cutoff probabilities called $\alpha$ and -$\beta$. +In {doc}`A Bayesian Formulation of Friedman and Wald's Problem ` +we presented a Bayesian procedure that makes +decisions by comparing a Bayesian posterior probability +$\pi$ with cutoff probabilities called $A$ and +$B$. To proceed, we borrow some Python code from the quantecon -lecture {doc}`A Problem that Stumped Milton Friedman ` -that computes $\alpha$ and $\beta$. +lecture {doc}`A Bayesian Formulation of Friedman and Wald's Problem ` +that computes optimal values of $A$ and $B$. ```{code-cell} python3 @jit(parallel=True) @@ -579,8 +579,8 @@ ax.plot(wf.π_grid, np.amin(np.column_stack([h_star, cost_L0, cost_L1]),axis=1), lw=15, alpha=0.1, color='b', label='minimum cost') -ax.annotate(r"$\beta$", xy=(β + 0.01, 0.5), fontsize=14) -ax.annotate(r"$\alpha$", xy=(α + 0.01, 0.5), fontsize=14) +ax.annotate(r"$B$", xy=(β + 0.01, 0.5), fontsize=14) +ax.annotate(r"$A$", xy=(α + 0.01, 0.5), fontsize=14) plt.vlines(β, 0, β * wf.L0, linestyle="--") plt.vlines(α, 0, (1 - α) * wf.L1, linestyle="--") @@ -595,14 +595,14 @@ plt.show() The above figure portrays the value function plotted against the decision maker’s Bayesian posterior. -It also shows the probabilities $\alpha$ and $\beta$. +It also shows the cutoff probabilities $A$ and $B$. The Bayesian decision rule is: -- accept $H_0$ if $\pi \geq \alpha$ -- accept $H_1$ if $\pi \leq \beta$ +- accept $H_0$ if $\pi \geq A$ +- accept $H_1$ if $\pi \leq B$ - delay deciding and draw another $z$ if - $\beta \leq \pi \leq \alpha$ + $ B \leq \pi \leq A$ We can calculate two “objective” loss functions under this situation conditioning on knowing for sure that nature has selected $f_{0}$, @@ -612,9 +612,9 @@ in the first case, or $f_{1}$, in the second case. $$ V^{0}\left(\pi\right)=\begin{cases} - 0 & \text{if }\alpha\leq\pi,\\ - c+EV^{0}\left(\pi^{\prime}\right) & \text{if }\beta\leq\pi<\alpha,\\ - \bar L_{1} & \text{if }\pi<\beta. + 0 & \text{if} A \leq\pi,\\ + c+EV^{0}\left(\pi^{\prime}\right) & \text{if }B\leq\pi< A,\\ + \bar L_{1} & \text{if }\pi= α] = wf.L0 + V[wf.π_grid >= A] = wf.L0 V_old = np.empty_like(V) while True: V_old[:] = V[:] - V[(β <= wf.π_grid) & (wf.π_grid < α)] = 0 + V[(B <= wf.π_grid) & (wf.π_grid < A)] = 0 for i in prange(len(wf.π_grid)): π = wf.π_grid[i] - if π >= α or π < β: + if π >= A or π < B: continue for j in prange(len(z_arr)): @@ -684,10 +684,10 @@ V1 = V_q(wf, 1) plt.plot(wf.π_grid, V0, label='$V^0$') plt.plot(wf.π_grid, V1, label='$V^1$') -plt.vlines(β, 0, wf.L0, linestyle='--') -plt.text(β+0.01, wf.L0/2, 'β') -plt.vlines(α, 0, wf.L0, linestyle='--') -plt.text(α+0.01, wf.L0/2, 'α') +plt.vlines(B, 0, wf.L0, linestyle='--') +plt.text(B+0.01, wf.L0/2, 'B') +plt.vlines(A, 0, wf.L0, linestyle='--') +plt.text(A+0.01, wf.L0/2, 'A') plt.xlabel(r'$\pi$') plt.title(r'Objective value function $V(\pi)$') plt.legend() @@ -778,8 +778,8 @@ plt.show() ## Was the Navy Captain’s Hunch Correct? -We now compare average (i.e., frequentist) losses obtained by the -frequentist and Bayesian decision rules. +We now compare average losses obtained by our frequentist Neyman-Pearson + and Bayesian decision rules. As a starting point, let’s compare average loss functions when $\pi^{*}=0.5$. @@ -806,7 +806,7 @@ plt.legend() plt.show() ``` -Evidently, there is no sample size $t$ at which the frequentist +Evidently, there is no sample size $t$ at which the Neyman-Pearson decision rule attains a lower loss function than does the Bayesian rule. Furthermore, the following graph indicates that the Bayesian decision @@ -841,7 +841,7 @@ $\pi^{*}=0.5=\pi_{0}$. π_star = 0.5 ``` -Recall that when $\pi^*=0.5$, the frequentist decision rule sets a +Recall that when $\pi^*=0.5$, the frequentist Neyman-Pearson decision rule sets a sample size `t_optimal` **ex ante**. For our parameter settings, we can compute its value: @@ -923,11 +923,11 @@ plt.show() ``` Later we’ll figure out how these distributions ultimately affect -objective expected values under the two decision rules. +objective expected values under the Neyman-Pearson and Bayesian decision rules. To begin, let’s look at simulations of the Bayesian’s beliefs over time. -We can easily compute the updated beliefs at any time $t$ using +We can compute updated beliefs at any time $t$ using the one-to-one mapping from $L_{t}$ to $\pi_{t}$ given $\pi_0$ described in this lecture {doc}`Likelihood Ratio Processes `. @@ -994,10 +994,10 @@ plt.show() ## Probability of Making Correct Decision -Now we use simulations to compute the fraction of samples in which the -Bayesian and the frequentist decision rules decide correctly. +Now we use simulations to compute the fractions of samples in which the +Bayesian and the frequentist Neyman-Pearson decision rules decide correctly. -For the frequentist rule, the probability of making the correct decision +For the frequentist Neyman-Pearson rule, the probability of making the correct decision under $f_{1}$ is the optimal probability of detection given $t$ that we defined earlier, and similarly it equals $1$ minus the optimal probability of a false alarm under $f_{0}$. @@ -1051,13 +1051,13 @@ plt.title('Uncond. probability of making correct decisions before t') plt.show() ``` -## Distribution of Likelihood Ratios at Frequentist’s $t$ +## Distribution of Likelihood Ratios at Neyman-Pearson's $t$ Next we use simulations to construct distributions of likelihood ratios after $t$ draws. To serve as useful reference points, we also show likelihood ratios that -correspond to the Bayesian cutoffs $\alpha$ and $\beta$. +correspond to the Bayesian cutoffs $A$ and $B$. In order to exhibit the distribution more clearly, we report logarithms of likelihood ratios. From 0050344a7f2ff70a54e01380dbf46738c4582901 Mon Sep 17 00:00:00 2001 From: Humphrey Yang Date: Wed, 6 Aug 2025 12:57:11 +1000 Subject: [PATCH 6/6] typo and code fixes --- lectures/navy_captain.md | 46 ++++++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/lectures/navy_captain.md b/lectures/navy_captain.md index 0846428b6..a7df21c53 100644 --- a/lectures/navy_captain.md +++ b/lectures/navy_captain.md @@ -327,7 +327,7 @@ plt.title("Receiver Operating Characteristic Curve") plt.show() ``` -We can minimizes the expected total loss presented in equation +We can minimize the expected total loss presented in equation {eq}`val1` by choosing $\left(t,d\right)$. Doing that delivers an expected loss @@ -355,7 +355,7 @@ def V_fre_d_t(d, t, L0_arr, L1_arr, π_star, wf): PFA = np.sum(L0_arr[:, t-1] < d) / N PD = np.sum(L1_arr[:, t-1] < d) / N - V = π_star * PFA *wf. L1 + (1 - π_star) * (1 - PD) * wf.L0 + V = π_star * PFA * wf.L1 + (1 - π_star) * (1 - PD) * wf.L0 return V ``` @@ -555,18 +555,18 @@ def find_cutoff_rule(wf, h): # The cutoff points can be found by differencing these costs with # The Bellman equation (J is always less than or equal to p_c_i) - β = π_grid[np.searchsorted( + B = π_grid[np.searchsorted( payoff_f1 - np.minimum(h, payoff_f0), 1e-10) - 1] - α = π_grid[np.searchsorted( + A = π_grid[np.searchsorted( np.minimum(h, payoff_f1) - payoff_f0, 1e-10) - 1] - return (β, α) + return (B, A) -β, α = find_cutoff_rule(wf, h_star) +B, A = find_cutoff_rule(wf, h_star) cost_L0 = (1 - wf.π_grid) * wf.L0 cost_L1 = wf.π_grid * wf.L1 @@ -579,11 +579,11 @@ ax.plot(wf.π_grid, np.amin(np.column_stack([h_star, cost_L0, cost_L1]),axis=1), lw=15, alpha=0.1, color='b', label='minimum cost') -ax.annotate(r"$B$", xy=(β + 0.01, 0.5), fontsize=14) -ax.annotate(r"$A$", xy=(α + 0.01, 0.5), fontsize=14) +ax.annotate(r"$B$", xy=(B + 0.01, 0.5), fontsize=14) +ax.annotate(r"$A$", xy=(A + 0.01, 0.5), fontsize=14) -plt.vlines(β, 0, β * wf.L0, linestyle="--") -plt.vlines(α, 0, (1 - α) * wf.L1, linestyle="--") +plt.vlines(B, 0, B * wf.L0, linestyle="--") +plt.vlines(A, 0, (1 - A) * wf.L1, linestyle="--") ax.set(xlim=(0, 1), ylim=(0, 0.5 * max(wf.L0, wf.L1)), ylabel="cost", xlabel=r"$\pi$", title="Value function") @@ -793,7 +793,7 @@ $\pi^{*}=0.5$. V_fre_arr, PFA_arr, PD_arr = compute_V_fre(L0_arr, L1_arr, π_star, wf) # bayesian -V_baye = π_star * V0 + π_star * V1 +V_baye = π_star * V0 + (1 - π_star) * V1 V_baye_bar = V_baye.min() ``` @@ -877,7 +877,7 @@ rule when $q= f_0$ and **later** when $q = f_1$. ```{code-cell} python3 @jit(parallel=True) -def check_results(L_arr, α, β, flag, π0): +def check_results(L_arr, A, B, flag, π0): N, T = L_arr.shape @@ -888,17 +888,17 @@ def check_results(L_arr, α, β, flag, π0): for i in prange(N): for t in range(T): - if (π_arr[i, t] < β) or (π_arr[i, t] > α): + if (π_arr[i, t] < B) or (π_arr[i, t] > A): time_arr[i] = t + 1 - correctness[i] = (flag == 0 and π_arr[i, t] > α) or (flag == 1 and π_arr[i, t] < β) + correctness[i] = (flag == 0 and π_arr[i, t] > A) or (flag == 1 and π_arr[i, t] < B) break return time_arr, correctness ``` ```{code-cell} python3 -time_arr0, correctness0 = check_results(L0_arr, α, β, 0, π_star) -time_arr1, correctness1 = check_results(L1_arr, α, β, 1, π_star) +time_arr0, correctness0 = check_results(L0_arr, A, B, 0, π_star) +time_arr1, correctness1 = check_results(L1_arr, A, B, 1, π_star) # unconditional distribution time_arr_u = np.concatenate((time_arr0, time_arr1)) @@ -963,7 +963,7 @@ The left graph compares $E\left(\pi_{t}\right)$ under $f_{0}$ to $1-E\left(\pi_{t}\right)$ under $f_{1}$: they lie on top of each other. -However, as the right hand size graph shows, there is significant +However, as the right hand side graph shows, there is significant difference in variances when $t$ is small: the variance is lower under $f_{1}$. @@ -1067,8 +1067,8 @@ $f_0$ generating the data, the other conditional on $f_1$ generating the data. ```{code-cell} python3 -Lα = (1 - π_star) * α / (π_star - π_star * α) -Lβ = (1 - π_star) * β / (π_star - π_star * β) +LA = (1 - π_star) * A / (π_star - π_star * A) +LB = (1 - π_star) * B / (π_star - π_star * B) ``` ```{code-cell} python3 @@ -1078,8 +1078,8 @@ bin_range = np.linspace(np.log(L_min), np.log(L_max), 50) n0 = plt.hist(np.log(L0_arr[:, t_idx]), bins=bin_range, alpha=0.4, label='f0 generates')[0] n1 = plt.hist(np.log(L1_arr[:, t_idx]), bins=bin_range, alpha=0.4, label='f1 generates')[0] -plt.vlines(np.log(Lβ), 0, max(n0.max(), n1.max()), linestyle='--', color='r', label='log($L_β$)') -plt.vlines(np.log(Lα), 0, max(n0.max(), n1.max()), linestyle='--', color='b', label='log($L_α$)') +plt.vlines(np.log(LB), 0, max(n0.max(), n1.max()), linestyle='--', color='r', label='log($L_B$)') +plt.vlines(np.log(LA), 0, max(n0.max(), n1.max()), linestyle='--', color='b', label='log($L_A$)') plt.legend() plt.xlabel('log(L)') @@ -1096,8 +1096,8 @@ distributions. ```{code-cell} python3 plt.hist(np.log(np.concatenate([L0_arr[:, t_idx], L1_arr[:, t_idx]])), bins=50, alpha=0.4, label='unconditional dist. of log(L)') -plt.vlines(np.log(Lβ), 0, max(n0.max(), n1.max()), linestyle='--', color='r', label='log($L_β$)') -plt.vlines(np.log(Lα), 0, max(n0.max(), n1.max()), linestyle='--', color='b', label='log($L_α$)') +plt.vlines(np.log(LB), 0, max(n0.max(), n1.max()), linestyle='--', color='r', label='log($L_B$)') +plt.vlines(np.log(LA), 0, max(n0.max(), n1.max()), linestyle='--', color='b', label='log($L_A$)') plt.legend() plt.xlabel('log(L)')