""" pu_auc.py ========= Reference implementation accompanying "Reliable AUC Evaluation for Positive-Unlabeled Classifiers: Calibrated Confidence Intervals under an Unknown Class Prior". Given a scoring model (anything with decision_function or predict_proba, in the scikit-learn convention) or precomputed real-valued scores, the routine returns 1. the decontaminated AUC against true negatives (point estimate), 2. a calibrated two-sided confidence interval (when the regular regime holds), or a one-sided lower bound near the singular boundary, 3. a degeneracy index I = n_U * d_hat^2 and a flag. Branches (the routine reports which it used and the guarantee status): bootstrap default robust branch. Resamples the whole pipeline, so the reported variance includes the uncertainty of the proportion AND of the score standardization. Distribution-free; applies to any scorer. Calibrated on the bi-Gaussian model in the paper. closed-form fast analytic branch (equations (5)-(11)). Valid when the achieved separation d_star is known; pass it through d_star, or accept that estimating it from a small test sample can degrade coverage. The plug-in variance is biased upward by the convex decontamination factor 1/(1-pi)^2 at high proportion and small unlabeled sample; the studentize option removes that bias. one-sided near the boundary the two-sided interval is replaced by a one-sided lower bound. Dependencies: numpy, scipy. No scikit-learn import is required. """ from __future__ import annotations from dataclasses import dataclass, field from typing import Callable, Optional, Sequence, Tuple, Union import numpy as np from scipy.stats import norm, multivariate_normal @dataclass class PUAUCResult: theta_hat: float theta_bc: float ci: Optional[Tuple[float, float]] lower_bound: Optional[float] pi_hat: float d_hat: float info_index: float regular: float flag: Union[int, str] method: str notes: str = "" components: dict = field(default_factory=dict) def __repr__(self) -> str: ci = "None" if self.ci is None else f"({self.ci[0]:.4f}, {self.ci[1]:.4f})" lb = "None" if self.lower_bound is None else f"{self.lower_bound:.4f}" return (f"PUAUCResult(theta_hat={self.theta_hat:.4f}, theta_bc={self.theta_bc:.4f}, " f"ci={ci}, lower_bound={lb}, pi_hat={self.pi_hat:.4f}, d_hat={self.d_hat:.4f}, " f"I={self.info_index:.2f}, flag={self.flag!r}, method={self.method!r})") # --------------------------------------------------------------------------- # # scoring helpers # # --------------------------------------------------------------------------- # def _score(model, X) -> np.ndarray: if hasattr(model, "decision_function"): return np.asarray(model.decision_function(X), dtype=float).ravel() if hasattr(model, "predict_proba"): proba = np.asarray(model.predict_proba(X), dtype=float) col = proba[:, -1] if proba.ndim == 2 else proba.ravel() col = np.clip(col, 1e-12, 1 - 1e-12) return np.log(col / (1.0 - col)) raise AttributeError("model must expose decision_function or predict_proba; " "or pass precomputed scores_pos and scores_unlab.") def _average_ranks(a: np.ndarray) -> np.ndarray: order = np.argsort(a, kind="mergesort") ranks = np.empty(a.size, dtype=float) s = a[order] i, n = 0, a.size while i < n: j = i while j + 1 < n and s[j + 1] == s[i]: j += 1 ranks[order[i:j + 1]] = 0.5 * (i + j) + 1.0 i = j + 1 return ranks def mann_whitney_auc(scores_pos: np.ndarray, scores_unlab: np.ndarray) -> float: sp = np.asarray(scores_pos, dtype=float) su = np.asarray(scores_unlab, dtype=float) n_p, n_u = sp.size, su.size ranks = _average_ranks(np.concatenate([sp, su])) return (ranks[:n_p].sum() - n_p * (n_p + 1) / 2.0) / (n_p * n_u) # --------------------------------------------------------------------------- # # proportion estimator # # --------------------------------------------------------------------------- # def _moment_mpe(scores_pos, scores_unlab, eps=1e-3, d_fixed=None): """ Two-moment estimator for a mixture with one known component. If d_fixed is given, the standardized separation is taken as known (the deterministic equivalent) and only the proportion is estimated; otherwise both are matched from the first two moments of the unlabeled sample. """ sp = np.asarray(scores_pos, dtype=float) su = np.asarray(scores_unlab, dtype=float) mu_p = sp.mean() sigma = sp.std(ddof=1) sigma = sigma if sigma > 0 else 1.0 n_u = su.size if d_fixed is None: m1, v1 = su.mean(), su.var(ddof=1) A = mu_p - m1 B = v1 - sigma ** 2 denom = A ** 2 + B pi0 = float(np.clip(B / denom if denom > 0 else eps, eps, 1 - eps)) g = A / (1.0 - pi0) d_hat = max(abs(g / sigma), 1e-6) mu_neg = mu_p - g else: d_hat = max(abs(float(d_fixed)), 1e-6) mu_neg = mu_p - d_hat * sigma z_pos = (sp - mu_neg) / sigma z_unlab = (su - mu_neg) / sigma pi_hat = float(np.clip(z_unlab.mean() / d_hat, eps, 1 - eps)) var_pi = (1.0 + pi_hat * (1.0 - pi_hat) * d_hat ** 2) / (n_u * d_hat ** 2) info = {"z_pos": z_pos, "z_unlab": z_unlab, "d_hat": d_hat, "mu_neg": mu_neg, "sigma": sigma} return pi_hat, var_pi, info def var_pi_two_channel(d_star, pi, n_p, n_u): """ Variance of the moment proportion estimator with the standardization estimated from the positive sample (Proposition 2). Two channels, Var(pi_hat) = a / n_U + b / n_P, a = (1 + pi (1 - pi) d^2) / d^2, # test channel, formula (9) b = 1 / d^2 + (1 - pi) ** 2 / 2. # standardization channel The test channel alone is formula (9), which conditions on the standardization being known. The standardization channel does not vanish when the positive sample is fixed and the unlabeled sample grows, and is often the larger of the two in deployment. Verified against Monte Carlo to within a few percent in the accompanying tests. """ d = max(abs(float(d_star)), 1e-6) a = (1.0 + pi * (1.0 - pi) * d ** 2) / d ** 2 b = 1.0 / d ** 2 + (1.0 - pi) ** 2 / 2.0 return a / n_u + b / n_p, a / n_u, b / n_p def _as_mpe(mpe, d_star): if mpe == "moments": return lambda sp, su: _moment_mpe(sp, su, d_fixed=d_star) if callable(mpe): return mpe raise ValueError("mpe must be 'moments' or a callable.") # --------------------------------------------------------------------------- # # closed-form bi-Gaussian variance # # --------------------------------------------------------------------------- # def _Phi2(x, y, rho): return float(multivariate_normal(mean=[0.0, 0.0], cov=[[1.0, rho], [rho, 1.0]]).cdf([x, y])) def _Psi(d): return _Phi2(d / np.sqrt(2.0), d / np.sqrt(2.0), 0.5) def _Psi0(d): return _Phi2(0.0, d / np.sqrt(2.0), 0.5) def _closed_form_var_dA(pi, d, n_p, n_u): theta_star = norm.cdf(d / np.sqrt(2.0)) mu_A = 0.5 * pi + (1.0 - pi) * theta_star sig_P2 = (pi ** 2) / 3.0 + 2.0 * pi * (1.0 - pi) * _Psi0(d) \ + (1.0 - pi) ** 2 * _Psi(d) - mu_A ** 2 sig_U2 = pi / 3.0 + (1.0 - pi) * _Psi(d) - mu_A ** 2 return max(sig_P2, 0.0) / n_p + max(sig_U2, 0.0) / n_u def _coupling_quadrature(pi, d, n_u, n_nodes=80): nodes, weights = np.polynomial.hermite_e.hermegauss(n_nodes) weights = weights / np.sqrt(2.0 * np.pi) def comp(mean): zz = nodes + mean h = norm.cdf(d - zz) psi = (zz - pi * d) / d return (np.sum(weights * h), np.sum(weights * psi), np.sum(weights * h * psi)) Eh_p, Ep_p, Ehp_p = comp(d) Eh_n, Ep_n, Ehp_n = comp(0.0) Eh = pi * Eh_p + (1 - pi) * Eh_n Ep = pi * Ep_p + (1 - pi) * Ep_n Ehp = pi * Ehp_p + (1 - pi) * Ehp_n return (Ehp - Eh * Ep) / n_u def _closed_form_V(pi_hat, d_hat, theta_hat, var_pi, n_p, n_u): var_dA = _closed_form_var_dA(pi_hat, d_hat, n_p, n_u) cov_q = _coupling_quadrature(pi_hat, d_hat, n_u) V = (var_dA + (theta_hat - 0.5) ** 2 * var_pi + 2.0 * (theta_hat - 0.5) * cov_q) / (1.0 - pi_hat) ** 2 return max(V, 0.0), var_dA, cov_q def _studentize_factor(var_pi, pi_hat): """ Second-order upward bias of the plug-in variance from the convex factor 1/(1-pi)^2: E[1/(1-pihat)^2] = (1/(1-pi)^2)(1 + 3 Var(pihat)/(1-pi)^2). Dividing the variance by this factor removes the leading bias. The factor is ~1 (no effect) when the proportion is well estimated, and grows at high proportion and small unlabeled sample, where it tightens the interval. """ return 1.0 + 3.0 * var_pi / (1.0 - pi_hat) ** 2 # --------------------------------------------------------------------------- # # bootstrap of the whole pipeline # # --------------------------------------------------------------------------- # def _bootstrap(scores_pos, scores_unlab, mpe_fn, n_boot, rng): sp = np.asarray(scores_pos, dtype=float) su = np.asarray(scores_unlab, dtype=float) n_p, n_u = sp.size, su.size thetas = np.empty(n_boot) pis = np.empty(n_boot) for b in range(n_boot): bp = sp[rng.integers(0, n_p, n_p)] bu = su[rng.integers(0, n_u, n_u)] pi_b, _, _ = mpe_fn(bp, bu) A_b = mann_whitney_auc(bp, bu) pis[b] = pi_b thetas[b] = np.clip((A_b - pi_b / 2.0) / (1.0 - pi_b), 0.0, 1.0) return thetas, pis # --------------------------------------------------------------------------- # # rigorous one-sided floor and Patra-Sen proportion bound # # --------------------------------------------------------------------------- # def _auc_with_sd(scores_pos, scores_unlab): """ Contaminated AUC mu_A = AUC(P,U) and its distribution-free DeLong standard deviation from empirical placement values. No model assumption. """ sp = np.asarray(scores_pos, dtype=float) su = np.asarray(scores_unlab, dtype=float) n_p, n_u = sp.size, su.size su_s = np.sort(su) sp_s = np.sort(sp) # placement of each positive against U, and each unlabeled against P, # with the mid-rank tie correction hP = (np.searchsorted(su_s, sp, side="right") + np.searchsorted(su_s, sp, side="left")) / (2.0 * n_u) hU = 1.0 - (np.searchsorted(sp_s, su, side="right") + np.searchsorted(sp_s, su, side="left")) / (2.0 * n_p) A = float(hP.mean()) sd = float(np.sqrt(hP.var(ddof=1) / n_p + hU.var(ddof=1) / n_u)) return A, sd def theta_floor_lower_bound(scores_pos, scores_unlab, alpha=0.05): """ Rigorous distribution-free one-sided lower bound on the decontaminated AUC. Since theta - mu_A = pi (theta - 1/2) >= 0 for theta >= 1/2, the target theta is at least the contaminated AUC mu_A, so a one-sided (1 - alpha) lower confidence limit for mu_A is a valid lower bound for theta. This needs no proportion estimate and no identifiability, which is why it is the honest report near the boundary where the proportion is not regularly estimable. """ A, sd = _auc_with_sd(scores_pos, scores_unlab) z1 = norm.ppf(1.0 - alpha) return float(np.clip(A - z1 * sd, 0.0, 1.0)), A, sd def patra_sen_lower_bound(scores_pos, scores_unlab, alpha_conf=0.05, gridsize=200): """ Honest distribution-free lower confidence bound on the weight of the unknown component a = 1 - pi (the negatives), in the spirit of Patra and Sen (2016). The known component is the positive distribution, estimated by the ECDF of the labeled-positive sample. For a candidate weight a, the implied unknown CDF is projected onto the set of valid CDFs (isotonic, clipped to [0,1]); the fitted mixture is compared to the unlabeled ECDF in L2, and the bound is the smallest a whose distance stays under a DKW threshold. The threshold is inflated for the estimation error of the known component, so the bound is conservative: P(a_lower <= 1 - pi) >= 1 - alpha_conf in the simulations of the accompanying tests. Near the boundary the bound degrades to near zero, which honestly reflects that few negatives can be certified there. Returns (a_lower, pi_upper) with pi_upper = 1 - a_lower, an upper confidence bound on the proportion (hence, through the identity, an upper bound on the decontaminated AUC under the model). """ su = np.asarray(scores_unlab, dtype=float) sp = np.asarray(scores_pos, dtype=float) n_u, n_p = su.size, sp.size su_s = np.sort(su) idx = np.linspace(0, n_u - 1, min(gridsize, n_u)).astype(int) grid = su_s[idx] Fn = np.searchsorted(su_s, grid, side="right") / n_u sp_s = np.sort(sp) Fb = np.searchsorted(sp_s, grid, side="right") / n_p c_beta = np.sqrt(np.log(2.0 / alpha_conf) / 2.0) thr = c_beta * (1.0 / np.sqrt(n_u) + 1.0 / np.sqrt(n_p)) def dist(a): Fs = (Fn - (1.0 - a) * Fb) / a proj = np.clip(_pava(Fs), 0.0, 1.0) fit = (1.0 - a) * Fb + a * proj return np.sqrt(np.mean((Fn - fit) ** 2)) a_lower = 1.0 for a in np.linspace(1.0, 2e-2, gridsize): if dist(a) <= thr: a_lower = float(a) return a_lower, float(1.0 - a_lower) def _pava(y): """Pool-adjacent-violators isotonic (nondecreasing) regression.""" y = np.asarray(y, dtype=float) vals, wts, lens = [], [], [] for i in range(y.size): v, w, ln = y[i], 1.0, 1 while vals and vals[-1] >= v: pv, pw, pl = vals.pop(), wts.pop(), lens.pop() v = (pw * pv + w * v) / (pw + w) w += pw ln += pl vals.append(v); wts.append(w); lens.append(ln) out = np.empty(y.size) pos = 0 for v, ln in zip(vals, lens): out[pos:pos + ln] = v pos += ln return out # --------------------------------------------------------------------------- # # main entry point # # --------------------------------------------------------------------------- # def pu_auc_inference( scores_pos: Optional[Sequence[float]] = None, scores_unlab: Optional[Sequence[float]] = None, *, model=None, X_pos=None, X_unlab=None, mpe: Union[str, Callable] = "moments", alpha: float = 0.05, method: str = "auto", d_star: Optional[float] = None, studentize: bool = True, n_boot: int = 500, info_low: float = 1.0, info_high: float = 4.0, width_warn: float = 0.5, sd_pi_warn: float = 0.10, sd_pi_degen: float = 0.25, random_state: Optional[int] = None, ) -> PUAUCResult: """ Parameters ---------- scores_pos, scores_unlab : 1-D score arrays, OR pass model with X_pos/X_unlab. mpe : "moments" or a callable mpe(sp, su) -> (pi_hat, var_pi, info) with info["d_hat"]. alpha : 1 - alpha is the target coverage. method : "auto"/"bootstrap" robust default, or "closed-form" fast analytic. d_star : if given, the standardized separation is treated as known (the deterministic equivalent), which makes the closed-form branch valid. studentize : remove the second-order bias of the closed-form plug-in variance from the convex decontamination factor (closed-form branch only). n_boot : bootstrap resamples for the robust branch. """ if scores_pos is None or scores_unlab is None: if model is None or X_pos is None or X_unlab is None: raise ValueError("Provide scores_pos/scores_unlab, or model with X_pos/X_unlab.") scores_pos = _score(model, X_pos) scores_unlab = _score(model, X_unlab) scores_pos = np.asarray(scores_pos, dtype=float).ravel() scores_unlab = np.asarray(scores_unlab, dtype=float).ravel() n_p, n_u = scores_pos.size, scores_unlab.size z = norm.ppf(1.0 - alpha / 2.0) rng = np.random.default_rng(random_state) mpe_fn = _as_mpe(mpe, d_star) # ---- full-sample fit -------------------------------------------------- # pi_hat, var_pi, info = mpe_fn(scores_pos, scores_unlab) d_hat = float(info["d_hat"]) A_obs = mann_whitney_auc(scores_pos, scores_unlab) theta_hat = float(np.clip((A_obs - pi_hat / 2.0) / (1.0 - pi_hat), 0.0, 1.0)) cov_quad = _coupling_quadrature(pi_hat, d_hat, n_u) bias = ((theta_hat - 0.5) * var_pi + cov_quad) / (1.0 - pi_hat) ** 2 theta_bc = float(np.clip(theta_hat - bias, 0.0, 1.0)) info_index = n_u * d_hat ** 2 # ---- closed-form fast branch ----------------------------------------- # if method == "closed-form": regular = float(np.clip( (np.log(max(info_index, 1e-9)) - np.log(info_low)) / (np.log(info_high) - np.log(info_low)), 0.0, 1.0)) V, var_dA, cov_q = _closed_form_V(pi_hat, d_hat, theta_hat, var_pi, n_p, n_u) stud = "" if studentize: V = V / _studentize_factor(var_pi, pi_hat) stud = (" Variance studentized: the second-order upward bias from the " "convex decontamination factor 1/(1-pi)^2 is removed.") if info_index < info_low: lb, A_obs2, sd_A = theta_floor_lower_bound(scores_pos, scores_unlab, alpha) a_low, pi_up = patra_sen_lower_bound(scores_pos, scores_unlab, alpha) return PUAUCResult(theta_hat, theta_bc, None, lb, pi_hat, d_hat, info_index, regular, 1, "one-sided", ("Degenerate regime (I below the crossover). The " "returned bound is the distribution-free floor " "theta >= mu_A; decontamination is not reliable " "here." + stud), {"V": V, "var_pi": var_pi, "mu_A": A_obs2, "sd_mu_A": sd_A, "patra_sen_a_lower": a_low, "pi_upper": pi_up}) half = z * np.sqrt(V) ci = (float(np.clip(theta_bc - half, 0.0, 1.0)), float(np.clip(theta_bc + half, 0.0, 1.0))) flag = 0 if info_index > info_high else "transition" note = ("Fast analytic branch (equations (5)-(11)). Valid when d_hat is " "known; pass d_star for the deterministic equivalent, otherwise " "estimating it from a small test sample can degrade coverage." + stud) return PUAUCResult(theta_hat, theta_bc, ci, None, pi_hat, d_hat, info_index, regular, flag, "closed-form", note, {"var_delta_A": var_dA, "var_pi": var_pi, "cov_delta_A_pi": cov_q, "V": V, "studentize_factor": _studentize_factor(var_pi, pi_hat)}) # ---- robust bootstrap default ---------------------------------------- # # The bootstrap dispersion of the proportion is the reliable degeneracy # signal; d_hat and I are not trustworthy near the boundary. thetas, pis = _bootstrap(scores_pos, scores_unlab, mpe_fn, n_boot, rng) V = float(thetas.var(ddof=1)) sd_pi = float(pis.std(ddof=1)) regular = float(np.clip(1.0 - (sd_pi - sd_pi_warn) / (sd_pi_degen - sd_pi_warn), 0.0, 1.0)) if sd_pi >= sd_pi_degen: lb, A_obs2, sd_A = theta_floor_lower_bound(scores_pos, scores_unlab, alpha) a_low, pi_up = patra_sen_lower_bound(scores_pos, scores_unlab, alpha) return PUAUCResult( theta_hat, theta_bc, None, lb, pi_hat, d_hat, info_index, regular, 1, "one-sided", (f"Degenerate regime: the proportion is weakly identified (bootstrap " f"sd of pi = {sd_pi:.2f}). A two-sided interval is not trustworthy. " "The returned bound is the distribution-free floor theta >= mu_A, " "which needs no proportion estimate; the Patra-Sen lower bound on " "1 - pi is reported as a diagnostic."), {"V": V, "boot_sd_pi": sd_pi, "n_boot": n_boot, "mu_A": A_obs2, "sd_mu_A": sd_A, "patra_sen_a_lower": a_low, "pi_upper": pi_up}) half = z * np.sqrt(V) ci = (float(np.clip(theta_bc - half, 0.0, 1.0)), float(np.clip(theta_bc + half, 0.0, 1.0))) note = ("Robust bootstrap branch. Interval is theta_bc +/- z * sd with sd the " "bootstrap standard deviation of the whole pipeline, so it captures " "the uncertainty of the proportion and of the score standardization. " "Distribution-free; calibrated on the bi-Gaussian model in the paper.") flag: Union[int, str] = 0 if sd_pi >= sd_pi_warn or (ci[1] - ci[0]) >= width_warn: flag = "transition" note += (f" Caution: the proportion is moderately identified (bootstrap " f"sd of pi = {sd_pi:.2f}); read the interval as approximate.") vp_two, vp_test, vp_std = var_pi_two_channel(d_hat, pi_hat, n_p, n_u) return PUAUCResult(theta_hat, theta_bc, ci, None, pi_hat, d_hat, info_index, regular, flag, "bootstrap", note, {"V": V, "boot_sd_theta": np.sqrt(V), "boot_sd_pi": sd_pi, "var_pi_formula": var_pi, "var_pi_two_channel": vp_two, "var_pi_test_channel": vp_test, "var_pi_std_channel": vp_std, "n_boot": n_boot}) # --------------------------------------------------------------------------- # # self-test # # --------------------------------------------------------------------------- # if __name__ == "__main__": rng = np.random.default_rng(0) # --- reliable branch: bootstrap coverage in the regular regime ---------- n_p, n_u = 400, 1600 d, pi = 1.74, 0.4 tt = norm.cdf(d / np.sqrt(2.0)) def draw(d, pi, n_u, rng): sp = rng.normal(d, 1.0, n_p) isp = rng.random(n_u) < pi su = np.where(isp, rng.normal(d, 1.0, n_u), rng.normal(0.0, 1.0, n_u)) return sp, su print("example:", pu_auc_inference(*draw(d, pi, n_u, rng), random_state=0)) cov = 0 reps = 80 for _ in range(reps): sp, su = draw(d, pi, n_u, rng) r = pu_auc_inference(sp, su, method="bootstrap", n_boot=300, random_state=1) cov += r.ci[0] <= tt <= r.ci[1] print(f"bootstrap coverage (reliable) d={d} pi={pi}: {cov / reps:.3f} (target 0.95)") # --- studentization validated in its proper regime ---------------------- # The studentize option removes the second-order bias of the closed-form # variance from the convex factor 1/(1-pi)^2. That effect is isolated when # the score standardization is known (negatives N(0,1), positives N(d,1)) # and only the proportion is estimated. We evaluate it directly on the # closed-form variance to avoid confounding with standardization estimation. def controlled_coverage(d, pi, n_p, n_u, studentize, reps, rng): zc = norm.ppf(0.975) tt = norm.cdf(d / np.sqrt(2.0)) hit = 0 for _ in range(reps): S_pos = rng.normal(d, 1.0, n_p) isp = rng.random(n_u) < pi S_U = np.where(isp, rng.normal(d, 1.0, n_u), rng.normal(0.0, 1.0, n_u)) A = mann_whitney_auc(S_pos, S_U) pih = float(np.clip(S_U.mean() / d, 1e-3, 1 - 1e-3)) # standardization known th = np.clip((A - pih / 2) / (1 - pih), 0.0, 1.0) varpi = (1 + pih * (1 - pih) * d ** 2) / (n_u * d ** 2) V, _, cq = _closed_form_V(pih, d, th, varpi, n_p, n_u) thbc = np.clip(th - ((th - 0.5) * varpi + cq) / (1 - pih) ** 2, 0.0, 1.0) if studentize: V = V / _studentize_factor(varpi, pih) hit += thbc - zc * np.sqrt(V) <= tt <= thbc + zc * np.sqrt(V) return hit / reps print("controlled regime (standardization known), closed-form coverage:") for cfg in [(1.74, 0.4, 400, 1600), (0.9, 0.6, 200, 400), (0.7, 0.7, 150, 300)]: d_, pi_, np_, nu_ = cfg cp = controlled_coverage(d_, pi_, np_, nu_, False, 4000, np.random.default_rng(7)) cs = controlled_coverage(d_, pi_, np_, nu_, True, 4000, np.random.default_rng(7)) print(f" d={d_} pi={pi_} n_U={nu_}: plain {cp:.3f} -> studentized {cs:.3f} (target 0.95)") # --- rigorous one-sided floor near the boundary ------------------------ print("one-sided floor theta >= mu_A near the boundary, P(theta >= bound):") for d_, pi_ in [(0.2, 0.5), (0.4, 0.4)]: tt = norm.cdf(d_ / np.sqrt(2.0)) rg = np.random.default_rng(11) hit = 0 reps = 1500 for _ in range(reps): sp = rg.normal(d_, 1.0, 400) isp = rg.random(1600) < pi_ su = np.where(isp, rg.normal(d_, 1.0, 1600), rg.normal(0.0, 1.0, 1600)) lb, _, _ = theta_floor_lower_bound(sp, su, alpha=0.05) hit += tt >= lb print(f" d={d_} pi={pi_}: {hit / reps:.3f} (target >= 0.95, rigorous floor)") # --- two-channel proportion variance (Proposition 2) ------------------- print("two-channel Var(pi_hat) = a/n_U + b/n_P vs Monte Carlo (std estimated):") for d_, pi_, np_, nu_ in [(1.74, 0.4, 400, 1600), (0.7, 0.5, 400, 1600)]: rg = np.random.default_rng(9) pis = np.empty(6000) for i in range(pis.size): sp = rg.normal(d_, 1.0, np_) isp = rg.random(nu_) < pi_ su = np.where(isp, rg.normal(d_, 1.0, nu_), rg.normal(0.0, 1.0, nu_)) pis[i], _, _ = _moment_mpe(sp, su, d_fixed=d_) pred, _, _ = var_pi_two_channel(d_, pi_, np_, nu_) print(f" d={d_} pi={pi_}: predicted {pred:.5f} MC {pis.var():.5f}")