@@ -946,44 +946,20 @@ $$ (eq:classerrorprob)
946946
947947where $\tilde \alpha_t = {\rm Prob}(l_t < 1 \mid f)$ and $\tilde \beta_t = {\rm Prob}(l_t \geq 1 \mid g)$.
948948
949- Now let's simulate timing protocol 2 and compute the classification error probability.
949+ Since for each $t$, the decision boundary is the same, the decision boundary can be computed as
950950
951951```{code-cell} ipython3
952- sequences_p2, true_sources_p2 = protocol_2(
953- π_minus_1, T_max, N_simulations)
954- l_ratios_p2, _ = compute_likelihood_ratios(sequences_p2)
955-
956- # Find decision boundary where f(w) = g(w)
957952root = brentq(lambda w: f(w) / g(w) - 1, 0.001, 0.999)
958-
959- # Compute theoretical tilde α_t and tilde β_t
960- def α_integrand(w):
961- """Integrand for tilde α_t = P(l_t < 1 | f)"""
962- return f(w) if f(w) / g(w) < 1 else 0
963-
964- def β_integrand(w):
965- """Integrand for tilde β_t = P(l_t >= 1 | g)"""
966- return g(w) if f(w) / g(w) >= 1 else 0
967-
968- # Compute the integrals
969- α_theory, _ = quad(α_integrand, 0, 1, limit=100)
970- β_theory, _ = quad(β_integrand, 0, 1, limit=100)
971-
972- theory_error = 0.5 * (α_theory + β_theory)
973-
974- print(f"theoretical tilde α_t = {α_theory:.4f}")
975- print(f"theoretical tilde β_t = {β_theory:.4f}")
976- print(f"theoretical classification error probability = {theory_error:.4f}")
977953```
978954
979- Since for each $t$, the decision boundary is the same, we can plot the distributions of $f$ and $g$ and the decision boundary
955+ we can plot the distributions of $f$ and $g$ and the decision boundary
980956
981957```{code-cell} ipython3
982958:tags: [hide-input]
983959
984960fig, ax = plt.subplots(figsize=(7, 6))
985961
986- w_range = np.linspace(0.001, 0.999 , 1000)
962+ w_range = np.linspace(1e-5, 1-1e-5 , 1000)
987963f_values = [f(w) for w in w_range]
988964g_values = [g(w) for w in w_range]
989965ratio_values = [f(w)/g(w) for w in w_range]
@@ -1025,13 +1001,40 @@ To the right of the green vertical line $g > f$, so $l_t >1 $; therefore a $
10251001
10261002 * The shaded blue area equals $\alpha$ -- the probability of classifying someone as a type $f$ when it is really a type $g$ individual.
10271003
1004+ This gives us clues about how to compute the theoretical classification error probability
10281005
1006+ ```{code-cell} ipython3
1007+ # Compute theoretical tilde α_t and tilde β_t
1008+ def α_integrand(w):
1009+ """Integrand for tilde α_t = P(l_t < 1 | f)"""
1010+ return f(w) if f(w) / g(w) < 1 else 0
10291011
1030- Let's see the classification algorithm performs in simulated data.
1012+ def β_integrand(w):
1013+ """Integrand for tilde β_t = P(l_t >= 1 | g)"""
1014+ return g(w) if f(w) / g(w) >= 1 else 0
1015+
1016+ # Compute the integrals
1017+ α_theory, _ = quad(α_integrand, 0, 1, limit=100)
1018+ β_theory, _ = quad(β_integrand, 0, 1, limit=100)
1019+
1020+ theory_error = 0.5 * (α_theory + β_theory)
1021+
1022+ print(f"theoretical tilde α_t = {α_theory:.4f}")
1023+ print(f"theoretical tilde β_t = {β_theory:.4f}")
1024+ print(f"theoretical classification error probability = {theory_error:.4f}")
1025+ ```
1026+
1027+ Now we simulate timing protocol 2 and compute the classification error probability.
1028+
1029+ In the next cell, we also compare the theoretical classification accuracy to the empirical classification accuracy
10311030
10321031```{code-cell} ipython3
10331032accuracy = np.empty(T_max)
10341033
1034+ sequences_p2, true_sources_p2 = protocol_2(
1035+ π_minus_1, T_max, N_simulations)
1036+ l_ratios_p2, _ = compute_likelihood_ratios(sequences_p2)
1037+
10351038for t in range(T_max):
10361039 predictions = (l_ratios_p2[:, t] >= 1)
10371040 actual = true_sources_p2[:, t]
@@ -1083,7 +1086,7 @@ We've already encountered one discrepancy measure -- the Kullback-Leibler (KL) d
10831086
10841087We now briefly explore two alternative discrepancy measures.
10851088
1086- ### Chernoff Entropy
1089+ ### Chernoff entropy
10871090
10881091Chernoff entropy was motivated by an early application of the [theory of large deviations](https://en.wikipedia.org/wiki/Large_deviations_theory).
10891092
10971100C(f,g) = - \log \min_ {\phi \in (0,1)} \int f^\phi(x) g^{1-\phi}(x) dx
10981101$$
10991102
1100- (TO TOM: since we are using natural log here, should we plot using $e^{-C(f,g)T}$ instead of $2^{-C(f,g)T}$?)
1101-
11021103An upper bound on model selection error probabilty is
11031104
11041105$$
11071108
11081109Thus, Chernoff entropy is an upper bound on the exponential rate at which the selection error probability falls as sample size $T$ grows.
11091110
1110- Let's compute Chernoff enropy numerically with some Python code
1111+ Let's compute Chernoff entropy numerically with some Python code
11111112
11121113```{code-cell} ipython3
11131114def chernoff_integrand(ϕ, f, g):
@@ -1117,7 +1118,7 @@ def chernoff_integrand(ϕ, f, g):
11171118 def integrand(w):
11181119 return f(w)**ϕ * g(w)**(1-ϕ)
11191120
1120- result, _ = quad(integrand, 0.0001, 0.9999 )
1121+ result, _ = quad(integrand, 1e-5, 1-1e-5 )
11211122 return result
11221123
11231124def compute_chernoff_entropy(f, g):
@@ -1130,7 +1131,7 @@ def compute_chernoff_entropy(f, g):
11301131 # Find the minimum over ϕ in (0,1)
11311132 result = minimize_scalar(objective,
11321133 # For numerical stability
1133- bounds=(0.001, 0.999 ),
1134+ bounds=(1e-5, 1-1e-5 ),
11341135 method='bounded')
11351136 min_value = result.fun
11361137 ϕ_optimal = result.x
@@ -1202,20 +1203,20 @@ def js_divergence(f, g):
12021203 def integrand(w):
12031204 ratio = p(w) / q(w)
12041205 return p(w) * np.log(ratio)
1205- result, _ = quad(integrand, 0.0001, 0.9999 )
1206+ result, _ = quad(integrand, 1e-5, 1-1e-5 )
12061207 return result
12071208
12081209 js_div = 0.5 * kl_div(f, m) + 0.5 * kl_div(g, m)
12091210 return js_div
12101211
12111212def kl_divergence(f, g):
12121213 """
1213- Compute KL divergence D (f, g)
1214+ Compute KL divergence KL (f, g)
12141215 """
12151216 def integrand(w):
12161217 return f(w) * np.log(f(w) / g(w))
12171218
1218- result, _ = quad(integrand, 0.0001, 0.9999 )
1219+ result, _ = quad(integrand, 1e-5, 1-1e-5 )
12191220 return result
12201221
12211222distribution_pairs = [
@@ -1384,7 +1385,7 @@ def error_divergence_cor():
13841385 N_sims = 5000
13851386 N_half = N_sims // 2
13861387
1387- # Initialize arrays to store results
1388+ # Initialize arrays
13881389 n_pairs = len(distribution_pairs)
13891390 kl_fg_vals = np.zeros(n_pairs)
13901391 kl_gf_vals = np.zeros(n_pairs)
@@ -1437,6 +1438,8 @@ cor_data = error_divergence_cor()
14371438Now let's visualize the correlations
14381439
14391440```{code-cell} ipython3
1441+ :tags: [hide-input]
1442+
14401443def plot_error_divergence(data):
14411444 """
14421445 Plot correlations between error probability and divergence measures.
@@ -1477,7 +1480,7 @@ plot_error_divergence(cor_data)
14771480
14781481Evidently, Chernoff entropy and Jensen-Shannon entropy each covary tightly with the model selection error probability.
14791482
1480- We'll see encounter related ideas in {doc}`wald_friedman`.
1483+ We'll encounter related ideas in {doc}`wald_friedman`.
14811484
14821485+++
14831486
0 commit comments