diff --git a/lectures/likelihood_ratio_process.md b/lectures/likelihood_ratio_process.md index a2037eae6..c9a7357bb 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,260 @@ 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_{ij}$ count 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)} +$$ + +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)}} +$$ (eq:llr_markov) + +### KL divergence rate + +By the ergodic theorem for irreducible, aperiodic Markov chains, we have + +$$ +\frac{N_{ij}}{T} \xrightarrow{a.s.} \pi_i^{(f)}P_{ij}^{(f)} \quad \text{as } T \to \infty +$$ + +where $\boldsymbol{\pi}^{(f)}$ is the stationary distribution satisfying $\boldsymbol{\pi}^{(f)} = \boldsymbol{\pi}^{(f)} P^{(f)}$. + +Therefore, + +$$ +\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)}} +$$ + +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)} \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)})} +$$ + +where $KL(P_{i\cdot}^{(f)}, P_{i\cdot}^{(g)})$ is the row-wise KL divergence. + + +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 +$$ + +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 +$$ + +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 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): + """ + 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 + """ + 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 + +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]]) + +# Compute stationary distributions +π_f = compute_stationary_dist(P_f) +π_g = compute_stationary_dist(P_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, π_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}") +``` + +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 line predicted by the theory + +```{code-cell} ipython3 +T = 500 +N_paths = 1000 +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, + π_f, π_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.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.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, + π_f, π_g) + +# Plot results +plt.figure(figsize=(10, 6)) +plt.plot(T_values, errors, linewidth=2) +plt.xlabel('$T$') +plt.ylabel('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 +2660,3 @@ $$ ```{solution-end} ``` - - - diff --git a/lectures/navy_captain.md b/lectures/navy_captain.md index 56431f556..a7df21c53 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 minimize the expected total loss presented in equation {eq}`val1` by choosing $\left(t,d\right)$. Doing that delivers an expected loss @@ -353,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 ``` @@ -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) @@ -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"$\beta$", xy=(β + 0.01, 0.5), fontsize=14) -ax.annotate(r"$\alpha$", 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") @@ -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$. @@ -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() ``` @@ -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: @@ -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)) @@ -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 `. @@ -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}$. @@ -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. @@ -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)')