Adding Bayesian Model (4 files)

bayesian.py

@@ -0,0 +1,1245 @@
+# bayesian.py
+# ============================================================
+# BAYESIAN FIT + CI + ELASTICITY
+# ============================================================
+
+import os
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+from scipy import optimize, stats
+from scipy.io import loadmat
+from scipy.optimize import brentq
+from scipy.special import betaln, gammaln
+
+from data_utils import get_data
+
+
+# ============================================================
+# 1) DATA LOADING
+# ============================================================
+
+def load_xy(
+    perc=95,
+    suv_path="suv_percentilesSLOthenUWM.mat",
+    flags_path="flags_combined.mat",
+):
+    """
+    Load feature x and binary label y.
+    """
+    here = os.path.dirname(os.path.abspath(__file__))
+
+    suv_full = os.path.join(here, suv_path)
+    flags_full = os.path.join(here, flags_path)
+
+    print("Loading SUV from:", suv_full)
+    print("Loading FLAGS from:", flags_full)
+
+    suv_dict = loadmat(suv_full)
+    flags_dict = loadmat(flags_full)
+
+    x, y = get_data(perc, suv_dict, flags_dict)
+
+    x = np.asarray(x, float).ravel()
+    y = np.asarray(y, int).ravel()
+
+    m = np.isfinite(x)
+    x, y = x[m], y[m]
+
+    x = np.clip(x, 1e-12, None)
+    return x, y
+
+
+# ============================================================
+# 2) CORE MODEL FUNCTIONS
+# ============================================================
+
+def logistic(z):
+    return 1.0 / (1.0 + np.exp(-np.clip(z, -60, 60)))
+
+
+def sigmoid(t):
+    return 1.0 / (1.0 + np.exp(-np.clip(t, -60, 60)))
+
+
+def dE_full(x, a, b, s, k, th):
+    """
+    log f_BP(x|a,b,s) - log f_Gamma(x|k,th), including constants.
+    """
+    x = np.asarray(x, float)
+    return (
+        (a - k) * np.log(x)
+        - (a + b) * np.log1p(x / s)
+        + x / th
+        - a * np.log(s)
+        - betaln(a, b)
+        + k * np.log(th)
+        + gammaln(k)
+    )
+
+
+def softplus(t):
+    t = np.asarray(t, float)
+    return np.log1p(np.exp(-np.abs(t))) + np.maximum(t, 0.0)
+
+
+def theta_max(a, b, k, s, eps=1e-12):
+    """
+    Monotonicity cap for theta.
+    """
+    A = a - k
+    if A <= 0:
+        return np.inf
+    r = np.sqrt(a + b) - np.sqrt(max(A, eps))
+    return np.inf if r <= 1e-12 else s / (r * r)
+
+
+def unpack(phi):
+    """
+    Reparameterisation:
+      phi = [p_raw, b_raw, s_raw, k_raw, d_raw, u_raw]
+
+      p in (0,1)
+      b,s,k > 0
+      a = k + delta with delta > 0
+      theta = theta_cap * sigmoid(u_raw)
+    """
+    p_raw, b_raw, s_raw, k_raw, d_raw, u_raw = phi
+
+    p = sigmoid(p_raw)
+
+    b = softplus(b_raw) + 1e-6
+    s = softplus(s_raw) + 1e-6
+    k = softplus(k_raw) + 1e-6
+
+    delta = softplus(d_raw) + 1e-6
+    a = k + delta
+
+    thcap = theta_max(a, b, k, s)
+    th = thcap * sigmoid(u_raw)
+
+    return p, a, b, s, k, th, thcap
+
+
+def make_priors(y, tau=25.0):
+    """
+    Beta(TAU*p_emp, TAU*(1-p_emp)) prior on prevalence p.
+    """
+    p_emp = float(np.mean(y))
+    alpha = max(tau * p_emp, 1e-6)
+    beta = max(tau * (1.0 - p_emp), 1e-6)
+    return alpha, beta
+
+
+def neg_post(phi, X, y, alpha, beta, use_prior_p=True, prior_r=(1.05, 1.05)):
+    """
+    Negative log-posterior = NLL + optional priors.
+
+    Priors used here:
+      - Beta prior on prevalence p
+      - Beta prior on r = theta/theta_cap
+
+    No extra priors on a, b, s, k.
+    """
+    p, a, b, s, k, th, thcap = unpack(phi)
+    eps = 1e-12
+
+    logit_val = (np.log(p) - np.log(1.0 - p)) + dE_full(X, a, b, s, k, th)
+    px = logistic(logit_val)
+
+    nll = -np.sum(y * np.log(px + eps) + (1 - y) * np.log(1 - px + eps))
+
+    if use_prior_p:
+        nll += -((alpha - 1) * np.log(p + eps) + (beta - 1) * np.log(1 - p + eps))
+
+    if prior_r is not None and np.isfinite(thcap) and thcap > 0:
+        r = np.clip(th / thcap, 1e-9, 1 - 1e-9)
+        nll += -((prior_r[0] - 1) * np.log(r) + (prior_r[1] - 1) * np.log(1 - r))
+
+    return float(nll)
+
+
+def init_phi(X, y):
+    """
+    Stable initial values.
+    """
+    X = np.asarray(X, float)
+    y = np.asarray(y, int)
+
+    X0 = X[y == 0]
+    m0 = X0.mean() if X0.size else X.mean()
+    v0 = X0.var() if X0.size else X.var()
+    k0 = 2.0 if v0 <= 0 else max((m0 * m0) / (v0 + 1e-9), 1.5)
+
+    X1 = X[y == 1]
+    m1 = np.median(X1) if X1.size else np.median(X)
+
+    p0 = np.clip(float(np.mean(y)), 1e-3, 1 - 1e-3)
+    b0 = 1.5
+    s0 = max(m1, 0.5)
+
+    return np.array(
+        [
+            np.log(p0 / (1 - p0)),         # p_raw
+            np.log(np.expm1(b0) + 1e-9),   # b_raw
+            np.log(np.expm1(s0) + 1e-9),   # s_raw
+            np.log(np.expm1(k0) + 1e-9),   # k_raw
+            np.log(np.expm1(1.0) + 1e-9),  # d_raw
+            -0.2,                          # u_raw
+        ],
+        dtype=float,
+    )
+
+
+def fit_bayes(X, y, seed=0, use_prior_p=True, prior_r=(1.05, 1.05), tau=25.0):
+    """
+    MAP fit using L-BFGS-B.
+    """
+    X = np.asarray(X, float)
+    y = np.asarray(y, int)
+
+    alpha, beta = make_priors(y, tau=tau)
+
+    obj = lambda w: neg_post(
+        w,
+        X,
+        y,
+        alpha,
+        beta,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+    )
+
+    w0 = init_phi(X, y)
+    res = optimize.minimize(
+        obj,
+        w0,
+        method="L-BFGS-B",
+        options={"maxiter": 6000, "ftol": 1e-9},
+    )
+
+    if not (res.success and np.isfinite(res.fun)):
+        rng = np.random.default_rng(seed)
+        w1 = w0 + rng.normal(0, 0.2, size=w0.shape)
+        res = optimize.minimize(
+            obj,
+            w1,
+            method="L-BFGS-B",
+            options={"maxiter": 6000, "ftol": 1e-9},
+        )
+
+    theta_hat = unpack(res.x)
+    return theta_hat, res
+
+
+def P_with(theta_hat, x):
+    """
+    Posterior risk curve P(AE|x) under fitted model.
+    """
+    p, a, b, s, k, th, _ = theta_hat
+    x = np.asarray(x, float)
+    logit_val = (np.log(p) - np.log(1 - p)) + dE_full(x, a, b, s, k, th)
+    return logistic(logit_val)
+
+
+# ============================================================
+# 3) ORIGINAL / TRIM DATASETS
+# ============================================================
+
+def make_trimmed_dataset(X, y, value_to_drop=2.48122597, tol=1e-3):
+    """
+    Remove point(s) with x approximately equal to value_to_drop.
+    """
+    X = np.asarray(X, float)
+    y = np.asarray(y, int)
+
+    mask_keep = np.abs(X - value_to_drop) > tol
+    removed_idx = np.where(~mask_keep)[0]
+
+    return {
+        "X_orig": X.copy(),
+        "y_orig": y.copy(),
+        "X_trim": X[mask_keep],
+        "y_trim": y[mask_keep],
+        "removed_idx": removed_idx,
+    }
+
+
+def run_bayesian_group_fit(
+    perc=95,
+    suv_path="suv_percentilesSLOthenUWM.mat",
+    flags_path="flags_combined.mat",
+    value_to_drop=2.48122597,
+    tol=1e-3,
+    use_prior_p=True,
+    prior_r=(1.05, 1.05),
+    tau=25.0,
+):
+    """
+    Load data, create ORIGINAL/TRIM datasets,
+    and fit constrained Bayesian group model on both.
+    """
+    X_all, y_all = load_xy(perc=perc, suv_path=suv_path, flags_path=flags_path)
+
+    ds = make_trimmed_dataset(X_all, y_all, value_to_drop=value_to_drop, tol=tol)
+
+    theta_orig, res_orig = fit_bayes(
+        ds["X_orig"],
+        ds["y_orig"],
+        seed=0,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+
+    theta_trim, res_trim = fit_bayes(
+        ds["X_trim"],
+        ds["y_trim"],
+        seed=1,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+
+    return {
+        **ds,
+        "theta_orig": theta_orig,
+        "theta_trim": theta_trim,
+        "res_orig": res_orig,
+        "res_trim": res_trim,
+    }
+
+
+def summarize_theta(theta_hat):
+    p, a, b, s, k, th, thcap = theta_hat
+    return {
+        "p": p,
+        "a": a,
+        "b": b,
+        "s": s,
+        "k": k,
+        "theta": th,
+        "theta_cap": thcap,
+        "r": th / thcap if np.isfinite(thcap) and thcap > 0 else np.nan,
+    }
+
+
+def plot_bayesian_orig_trim_raw(
+    fit_results,
+    xmax=6.0,
+    suptitle="Conditional Probability of AE",
+    figsize=(12, 7),
+    dpi=140,
+):
+    """
+    One raw-x plot:
+      - ORIGINAL curve
+      - TRIM curve
+      - ORIGINAL data dots
+      - highlight removed point(s)
+    """
+    X_orig = fit_results["X_orig"]
+    y_orig = fit_results["y_orig"]
+    removed_idx = fit_results["removed_idx"]
+    theta_orig = fit_results["theta_orig"]
+    theta_trim = fit_results["theta_trim"]
+
+    fig, ax = plt.subplots(figsize=figsize, dpi=dpi)
+    fig.text(0.02, 0.5, suptitle, va="center", rotation="vertical", fontsize=14)
+
+    x_min = max(float(np.min(X_orig)), 1e-8)
+    x_max = float(xmax)
+
+    ax.set_xlim(x_min, x_max)
+    ax.set_xlabel("x")
+    ax.set_ylabel("P(AE | x)")
+    ax.set_ylim(-0.10, 1.10)
+    ax.grid(alpha=0.35)
+
+    x_grid = np.exp(np.linspace(np.log(x_min), np.log(x_max), 900))
+    p_curve_orig = P_with(theta_orig, x_grid)
+    p_curve_trim = P_with(theta_trim, x_grid)
+
+    l1, = ax.plot(x_grid, p_curve_orig, lw=2.2, color="C0", label="ORIGINAL (Bayesian fit)")
+    l2, = ax.plot(x_grid, p_curve_trim, lw=2.2, color="C1", label="TRIM (Bayesian fit)")
+
+    rng = np.random.default_rng(999)
+    jit = (rng.random(len(y_orig)) - 0.5) * 0.06
+
+    d_nc = ax.scatter(
+        X_orig[y_orig == 0],
+        (y_orig + jit)[y_orig == 0],
+        s=22,
+        alpha=0.65,
+        edgecolors="none",
+        color="C0",
+        label="NC samples (ORIGINAL)",
+    )
+    d_ae = ax.scatter(
+        X_orig[y_orig == 1],
+        (y_orig + jit)[y_orig == 1],
+        s=26,
+        alpha=0.85,
+        edgecolors="none",
+        color="C1",
+        label="AE samples (ORIGINAL)",
+    )
+
+    dout = None
+    if removed_idx.size > 0:
+        for j, i in enumerate(removed_idx):
+            jit_out = (rng.random() - 0.5) * 0.06
+            label = "Removed point" if j == 0 else None
+            dout = ax.scatter(
+                [float(X_orig[i])],
+                [float(y_orig[i] + jit_out)],
+                marker="x",
+                s=90,
+                linewidths=2,
+                color="k",
+                label=label,
+            )
+
+    handles = [l1, l2, d_nc, d_ae]
+    if dout is not None:
+        handles.append(dout)
+    labels = [h.get_label() for h in handles]
+    ax.legend(handles, labels, frameon=False, ncol=2, loc="lower right")
+
+    plt.tight_layout(rect=(0.06, 0.0, 1.0, 1.0))
+    return fig, ax
+
+
+# ============================================================
+# 4) CI ESTIMATION
+# ============================================================
+
+def x_at_p(theta_hat, p_target=0.5, lo=1e-6, hi=10.0):
+    """
+    Solve P(AE|x) = p_target for x.
+    """
+    f = lambda x: P_with(theta_hat, x) - p_target
+    try:
+        if f(lo) * f(hi) > 0:
+            return np.nan
+        return float(brentq(f, lo, hi))
+    except Exception:
+        return np.nan
+
+
+def slope_at_x(theta_hat, x0):
+    """
+    Numerical derivative of P(AE|x) at x0.
+    """
+    if not np.isfinite(x0):
+        return np.nan
+    h = 1e-3 * (1 + abs(x0))
+    return float((P_with(theta_hat, x0 + h) - P_with(theta_hat, x0 - h)) / (2 * h))
+
+
+def hess_fd(F, x):
+    """
+    Finite-difference Hessian.
+    """
+    x = np.asarray(x, float)
+    n = x.size
+    H = np.zeros((n, n))
+    h = 1e-4 * (1 + np.abs(x))
+
+    def grad_fd(G, z):
+        g = np.zeros_like(z)
+        for j in range(n):
+            ej = np.zeros_like(z)
+            ej[j] = 1.0
+            g[j] = (G(z + h[j] * ej) - G(z - h[j] * ej)) / (2 * h[j])
+        return g
+
+    for i in range(n):
+        ei = np.zeros_like(x)
+        ei[i] = 1.0
+        g_plus = grad_fd(F, x + h[i] * ei)
+        g_minus = grad_fd(F, x - h[i] * ei)
+        H[:, i] = (g_plus - g_minus) / (2 * h[i])
+
+    return 0.5 * (H + H.T)
+
+
+def jac_fd(Fvec, w):
+    """
+    Finite-difference Jacobian for vector-valued function.
+    """
+    f0 = Fvec(w)
+    m = f0.size
+    n = w.size
+    J = np.zeros((m, n))
+    h = 1e-4 * (1 + np.abs(w))
+
+    for j in range(n):
+        ej = np.zeros_like(w)
+        ej[j] = 1.0
+        J[:, j] = (Fvec(w + h[j] * ej) - Fvec(w - h[j] * ej)) / (2 * h[j])
+
+    return J
+
+
+def estimate_ci_bundle(
+    X,
+    y,
+    label,
+    x_grid,
+    B_nonpar=400,
+    B_param=400,
+    seed=123,
+    use_prior_p=True,
+    prior_r=(1.05, 1.05),
+    tau=25.0,
+):
+    X = np.asarray(X, float)
+    y = np.asarray(y, int)
+    rng = np.random.default_rng(seed)
+    n = len(X)
+
+    theta_hat, res = fit_bayes(
+        X,
+        y,
+        seed=seed,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+
+    pmap = P_with(theta_hat, x_grid)
+    x50 = x_at_p(theta_hat, 0.5, lo=max(1e-6, x_grid.min()), hi=x_grid.max())
+    s50 = slope_at_x(theta_hat, x50)
+
+    phi_hat = res.x
+    alpha, beta = make_priors(y, tau=tau)
+
+    H = hess_fd(
+        lambda w: neg_post(
+            w,
+            X,
+            y,
+            alpha,
+            beta,
+            use_prior_p=use_prior_p,
+            prior_r=prior_r,
+        ),
+        phi_hat,
+    )
+
+    Jp = jac_fd(lambda w: P_with(unpack(w), x_grid), phi_hat)
+
+    try:
+        Sigma_phi = np.linalg.inv(H)
+    except np.linalg.LinAlgError:
+        Sigma_phi = np.linalg.pinv(H)
+
+    var_p = np.einsum("ij,jk,ik->i", Jp, Sigma_phi, Jp)
+    se_p = np.sqrt(np.maximum(var_p, 0.0))
+    wald_lo = np.clip(pmap - 1.96 * se_p, 0, 1)
+    wald_hi = np.clip(pmap + 1.96 * se_p, 0, 1)
+
+    def theta_vec_from_phi(w):
+        p, a, b, s, k, th, thcap = unpack(w)
+        r = th / thcap if np.isfinite(thcap) and thcap > 0 else np.nan
+        return np.array([p, a, b, s, k, th, r], float)
+
+    Jtheta = jac_fd(theta_vec_from_phi, phi_hat)
+    Sigma_theta = Jtheta @ Sigma_phi @ Jtheta.T
+    theta_hat_vec = theta_vec_from_phi(phi_hat)
+    se_theta = np.sqrt(np.maximum(np.diag(Sigma_theta), 0.0))
+    wald_param_lo = theta_hat_vec - 1.96 * se_theta
+    wald_param_hi = theta_hat_vec + 1.96 * se_theta
+
+    curves_np = []
+    theta_np = []
+    x50_np = []
+    used_np = 0
+
+    for _ in range(B_nonpar):
+        idx = rng.integers(0, n, n)
+        Xb, yb = X[idx], y[idx]
+
+        if yb.sum() == 0 or yb.sum() == len(yb):
+            continue
+
+        try:
+            thb, rb = fit_bayes(
+                Xb,
+                yb,
+                seed=int(rng.integers(0, 10_000_000)),
+                use_prior_p=use_prior_p,
+                prior_r=prior_r,
+                tau=tau,
+            )
+            if not rb.success or not np.isfinite(rb.fun):
+                continue
+
+            curves_np.append(P_with(thb, x_grid))
+
+            p_b, a_b, b_b, s_b, k_b, th_b, thcap_b = thb
+            r_b = th_b / thcap_b if np.isfinite(thcap_b) and thcap_b > 0 else np.nan
+            theta_np.append([p_b, a_b, b_b, s_b, k_b, th_b, r_b])
+
+            x50_b = x_at_p(thb, 0.5, lo=max(1e-6, x_grid.min()), hi=x_grid.max())
+            x50_np.append(x50_b)
+
+            used_np += 1
+        except Exception:
+            continue
+
+    curves_np = np.asarray(curves_np)
+    theta_np = np.asarray(theta_np, float) if len(theta_np) else np.empty((0, 7))
+    x50_np = np.asarray(x50_np, float) if len(x50_np) else np.empty((0,))
+
+    np_lo = np.percentile(curves_np, 2.5, axis=0) if used_np else None
+    np_hi = np.percentile(curves_np, 97.5, axis=0) if used_np else None
+
+    curves_pb = []
+    theta_pb = []
+    x50_pb = []
+    used_pb = 0
+
+    p_hat, a_hat, b_hat, s_hat, k_hat, th_hat, _ = theta_hat
+
+    for _ in range(B_param):
+        yb = rng.binomial(1, p_hat, size=n)
+
+        if yb.sum() == 0 or yb.sum() == n:
+            continue
+
+        Xb = np.zeros(n, dtype=float)
+
+        idx_nc = np.where(yb == 0)[0]
+        idx_ae = np.where(yb == 1)[0]
+
+        if len(idx_nc) > 0:
+            Xb[idx_nc] = stats.gamma.rvs(
+                k_hat,
+                scale=th_hat,
+                size=len(idx_nc),
+                random_state=rng,
+            )
+
+        if len(idx_ae) > 0:
+            Xb[idx_ae] = stats.betaprime.rvs(
+                a_hat,
+                b_hat,
+                scale=s_hat,
+                size=len(idx_ae),
+                random_state=rng,
+            )
+
+        Xb = np.clip(Xb, 1e-12, None)
+
+        try:
+            thb, rb = fit_bayes(
+                Xb,
+                yb,
+                seed=int(rng.integers(0, 10_000_000)),
+                use_prior_p=use_prior_p,
+                prior_r=prior_r,
+                tau=tau,
+            )
+            if not rb.success or not np.isfinite(rb.fun):
+                continue
+
+            curves_pb.append(P_with(thb, x_grid))
+
+            p_b, a_b, b_b, s_b, k_b, th_b, thcap_b = thb
+            r_b = th_b / thcap_b if np.isfinite(thcap_b) and thcap_b > 0 else np.nan
+            theta_pb.append([p_b, a_b, b_b, s_b, k_b, th_b, r_b])
+
+            x50_b = x_at_p(thb, 0.5, lo=max(1e-6, x_grid.min()), hi=x_grid.max())
+            x50_pb.append(x50_b)
+
+            used_pb += 1
+        except Exception:
+            continue
+
+    curves_pb = np.asarray(curves_pb)
+    theta_pb = np.asarray(theta_pb, float) if len(theta_pb) else np.empty((0, 7))
+    x50_pb = np.asarray(x50_pb, float) if len(x50_pb) else np.empty((0,))
+
+    pb_lo = np.percentile(curves_pb, 2.5, axis=0) if used_pb else None
+    pb_hi = np.percentile(curves_pb, 97.5, axis=0) if used_pb else None
+
+    return {
+        "label": label,
+        "theta_hat": theta_hat,
+        "res": res,
+        "x50": x50,
+        "s50": s50,
+        "pmap": pmap,
+        "wald_lo": wald_lo,
+        "wald_hi": wald_hi,
+        "np_lo": np_lo,
+        "np_hi": np_hi,
+        "pb_lo": pb_lo,
+        "pb_hi": pb_hi,
+        "used_np": used_np,
+        "used_pb": used_pb,
+        "theta_hat_vec": theta_hat_vec,
+        "wald_param_lo": wald_param_lo,
+        "wald_param_hi": wald_param_hi,
+        "theta_np": theta_np,
+        "theta_pb": theta_pb,
+        "x50_np": x50_np,
+        "x50_pb": x50_pb,
+    }
+
+
+def run_bayesian_ci(
+    perc=95,
+    suv_path="suv_percentilesSLOthenUWM.mat",
+    flags_path="flags_combined.mat",
+    value_to_drop=2.48122597,
+    tol=1e-3,
+    xmax=10.0,
+    n_grid=600,
+    B_nonpar=400,
+    B_param=400,
+    seed=123,
+    use_prior_p=True,
+    prior_r=(1.05, 1.05),
+    tau=25.0,
+):
+    """
+    Run CI estimation for both ORIGINAL and TRIM datasets.
+    """
+    X_all, y_all = load_xy(perc=perc, suv_path=suv_path, flags_path=flags_path)
+    ds = make_trimmed_dataset(X_all, y_all, value_to_drop=value_to_drop, tol=tol)
+
+    x_grid = np.linspace(0, xmax, n_grid)
+
+    out_orig = estimate_ci_bundle(
+        ds["X_orig"],
+        ds["y_orig"],
+        label="ORIGINAL",
+        x_grid=x_grid,
+        B_nonpar=B_nonpar,
+        B_param=B_param,
+        seed=seed,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+
+    out_trim = estimate_ci_bundle(
+        ds["X_trim"],
+        ds["y_trim"],
+        label="TRIM",
+        x_grid=x_grid,
+        B_nonpar=B_nonpar,
+        B_param=B_param,
+        seed=seed + 1,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+
+    return {
+        **ds,
+        "x_grid": x_grid,
+        "orig": out_orig,
+        "trim": out_trim,
+    }
+
+
+def make_param_ci_table(ci_res):
+    """
+    Parameter CI table for FULL and TRIM, for Wald / Nonparam / Parametric.
+    """
+    rows = []
+    names = ["p", "a", "b", "s", "k", "theta", "r"]
+
+    for dataset_key, dataset_name in [("orig", "FULL"), ("trim", "TRIM")]:
+        out = ci_res[dataset_key]
+        hat = out["theta_hat_vec"]
+
+        for i, name in enumerate(names):
+            rows.append(
+                {
+                    "Dataset": dataset_name,
+                    "Method": "Wald",
+                    "Parameter": name,
+                    "Estimate": hat[i],
+                    "LL": out["wald_param_lo"][i],
+                    "UL": out["wald_param_hi"][i],
+                }
+            )
+
+        if out["theta_np"].shape[0] > 0:
+            lo = np.nanpercentile(out["theta_np"], 2.5, axis=0)
+            hi = np.nanpercentile(out["theta_np"], 97.5, axis=0)
+            for i, name in enumerate(names):
+                rows.append(
+                    {
+                        "Dataset": dataset_name,
+                        "Method": "Nonparam",
+                        "Parameter": name,
+                        "Estimate": hat[i],
+                        "LL": lo[i],
+                        "UL": hi[i],
+                    }
+                )
+
+        if out["theta_pb"].shape[0] > 0:
+            lo = np.nanpercentile(out["theta_pb"], 2.5, axis=0)
+            hi = np.nanpercentile(out["theta_pb"], 97.5, axis=0)
+            for i, name in enumerate(names):
+                rows.append(
+                    {
+                        "Dataset": dataset_name,
+                        "Method": "Parametric",
+                        "Parameter": name,
+                        "Estimate": hat[i],
+                        "LL": lo[i],
+                        "UL": hi[i],
+                    }
+                )
+
+    return pd.DataFrame(rows)
+
+
+def make_x50_ci_table(ci_res):
+    """
+    x50 CI table for FULL and TRIM, for Wald / Nonparam / Parametric.
+    """
+    rows = []
+
+    for dataset_key, dataset_name in [("orig", "FULL"), ("trim", "TRIM")]:
+        out = ci_res[dataset_key]
+
+        wald_x50_lo = np.nan
+        wald_x50_hi = np.nan
+        try:
+            wald_x50_lo = np.interp(0.5, out["wald_lo"], ci_res["x_grid"])
+            wald_x50_hi = np.interp(0.5, out["wald_hi"], ci_res["x_grid"])
+        except Exception:
+            pass
+
+        rows.append(
+            {
+                "Dataset": dataset_name,
+                "Method": "Wald",
+                "Estimate": out["x50"],
+                "LL": wald_x50_lo,
+                "UL": wald_x50_hi,
+            }
+        )
+
+        if len(out["x50_np"]) > 0:
+            rows.append(
+                {
+                    "Dataset": dataset_name,
+                    "Method": "Nonparam",
+                    "Estimate": out["x50"],
+                    "LL": np.nanpercentile(out["x50_np"], 2.5),
+                    "UL": np.nanpercentile(out["x50_np"], 97.5),
+                }
+            )
+
+        if len(out["x50_pb"]) > 0:
+            rows.append(
+                {
+                    "Dataset": dataset_name,
+                    "Method": "Parametric",
+                    "Estimate": out["x50"],
+                    "LL": np.nanpercentile(out["x50_pb"], 2.5),
+                    "UL": np.nanpercentile(out["x50_pb"], 97.5),
+                }
+            )
+
+    return pd.DataFrame(rows)
+
+
+def make_curve_ci_table(ci_res, grid_every=25):
+    """
+    Long-format curve CI table.
+    Contains LL/UL of P(AE|X) across x-grid for all methods.
+    """
+    rows = []
+    x_grid = ci_res["x_grid"][::grid_every]
+
+    for dataset_key, dataset_name in [("orig", "FULL"), ("trim", "TRIM")]:
+        out = ci_res[dataset_key]
+
+        for method, lo_key, hi_key in [
+            ("Wald", "wald_lo", "wald_hi"),
+            ("Nonparam", "np_lo", "np_hi"),
+            ("Parametric", "pb_lo", "pb_hi"),
+        ]:
+            lo = out.get(lo_key, None)
+            hi = out.get(hi_key, None)
+            if lo is None or hi is None:
+                continue
+
+            lo = lo[::grid_every]
+            hi = hi[::grid_every]
+            est = out["pmap"][::grid_every]
+
+            for x, e, l, u in zip(x_grid, est, lo, hi):
+                rows.append(
+                    {
+                        "Dataset": dataset_name,
+                        "Method": method,
+                        "X": x,
+                        "Estimate": e,
+                        "LL": l,
+                        "UL": u,
+                    }
+                )
+
+    return pd.DataFrame(rows)
+
+
+def plot_ci_original_trim(ci_res, xmax=6.0):
+    """
+    2-panel figure:
+      left  = FULL
+      right = TRIM
+    """
+    x_grid = ci_res["x_grid"]
+    rng = np.random.default_rng(999)
+
+    fig, axes = plt.subplots(1, 2, figsize=(16, 7.0), dpi=150, sharey=True)
+
+    for ax, X, y, out, title, panel in [
+        (axes[0], ci_res["X_orig"], ci_res["y_orig"], ci_res["orig"], "FULL", "A"),
+        (axes[1], ci_res["X_trim"], ci_res["y_trim"], ci_res["trim"], "TRIM", "B"),
+    ]:
+        ax.fill_between(x_grid, out["wald_lo"], out["wald_hi"], color="#2ca02c", alpha=0.10)
+
+        if out["used_np"]:
+            ax.fill_between(x_grid, out["np_lo"], out["np_hi"], color="#17becf", alpha=0.10)
+
+        if out["used_pb"]:
+            ax.fill_between(x_grid, out["pb_lo"], out["pb_hi"], color="#e91e63", alpha=0.10)
+
+        ax.plot(x_grid, out["wald_lo"], color="#2ca02c", lw=1.6, ls="--")
+        h_wald, = ax.plot(
+            x_grid,
+            out["wald_hi"],
+            color="#2ca02c",
+            lw=1.6,
+            ls="--",
+            label="CI: Wald (delta) 95%",
+        )
+
+        h_np = None
+        if out["used_np"]:
+            ax.plot(x_grid, out["np_lo"], color="#17becf", lw=1.6, ls=(0, (1, 2)))
+            h_np, = ax.plot(
+                x_grid,
+                out["np_hi"],
+                color="#17becf",
+                lw=1.6,
+                ls=(0, (1, 2)),
+                label="CI: Nonparam bootstrap 95%",
+            )
+
+        h_pb = None
+        if out["used_pb"]:
+            ax.plot(x_grid, out["pb_lo"], color="#e91e63", lw=1.6, ls="-.")
+            h_pb, = ax.plot(
+                x_grid,
+                out["pb_hi"],
+                color="#e91e63",
+                lw=1.6,
+                ls="-.",
+                label="CI: Parametric bootstrap 95%",
+            )
+
+        h_fit, = ax.plot(x_grid, out["pmap"], color="k", lw=2.4, label="Bayesian fit")
+
+        jit = (rng.random(len(y)) - 0.5) * 0.035
+
+        h_nc = ax.scatter(
+            X[y == 0],
+            (y + jit)[y == 0],
+            s=18,
+            alpha=0.55,
+            color="#5dade2",
+            edgecolors="none",
+            label="NC",
+        )
+
+        h_ae = ax.scatter(
+            X[y == 1],
+            (y + jit)[y == 1],
+            s=24,
+            alpha=0.80,
+            color="#f39c3d",
+            edgecolors="none",
+            label="AE",
+        )
+
+        ax.text(
+            0.02,
+            0.98,
+            panel,
+            transform=ax.transAxes,
+            ha="left",
+            va="top",
+            fontsize=16,
+            fontweight="bold",
+        )
+
+        ax.set_xlim(0, xmax)
+        ax.set_ylim(-0.05, 1.05)
+        ax.set_xlabel("X", fontsize=13, fontweight="bold")
+        ax.set_title(title, fontsize=13, fontweight="bold")
+        ax.grid(alpha=0.25)
+
+        handles = [h_nc, h_ae, h_fit, h_wald]
+        if h_np is not None:
+            handles.append(h_np)
+        if h_pb is not None:
+            handles.append(h_pb)
+
+        labels = [h.get_label() for h in handles]
+
+        ax.legend(
+            handles,
+            labels,
+            loc="upper center",
+            bbox_to_anchor=(0.5, -0.20),
+            ncol=2,
+            frameon=False,
+            fontsize=10,
+        )
+
+    axes[0].set_ylabel("P(AE | X)", fontsize=13, fontweight="bold")
+    plt.tight_layout(rect=(0, 0.08, 1, 1))
+    return fig, axes
+
+
+# ============================================================
+# 5) X50 ELASTICITY ANALYSIS
+# ============================================================
+
+PARAM_NAMES_X50_ELAS = [r"$\pi$", r"$a$", r"$b$", r"$s$", r"$k$", r"$\vartheta$"]
+
+
+def theta6_from_hat(theta_hat):
+    """
+    Extract the first 6 raw model parameters from theta_hat:
+      (p, a, b, s, k, th)
+    """
+    th = np.asarray(theta_hat, float).ravel()
+    if th.size < 6:
+        raise ValueError(f"Expected at least 6 parameters, got {th.size}")
+    return th[:6].copy()
+
+
+def step_vec_theta(theta, rel_step=1e-6, abs_min=1e-10):
+    """
+    Relative finite-difference step on raw theta scale.
+    """
+    theta = np.asarray(theta, float).ravel()
+    return np.maximum(abs_min, rel_step * np.maximum(1.0, np.abs(theta)))
+
+
+def grad_central_theta(F_theta, theta0, rel_step=1e-6, abs_min=1e-10, pi_eps=1e-12):
+    """
+    Central differences in RAW theta.
+    Keeps:
+      - p in (pi_eps, 1-pi_eps)
+      - positive parameters > 0
+    """
+    theta0 = np.asarray(theta0, float).ravel()
+    h = step_vec_theta(theta0, rel_step=rel_step, abs_min=abs_min)
+    g = np.zeros_like(theta0)
+
+    for j in range(theta0.size):
+        th_plus = theta0.copy()
+        th_minus = theta0.copy()
+        hj = h[j]
+
+        if j == 0:
+            p0 = float(np.clip(theta0[0], pi_eps, 1 - pi_eps))
+            hj = min(hj, p0 - pi_eps, (1 - pi_eps) - p0)
+            hj = max(hj, abs_min)
+            th_plus[0] = np.clip(p0 + hj, pi_eps, 1 - pi_eps)
+            th_minus[0] = np.clip(p0 - hj, pi_eps, 1 - pi_eps)
+        else:
+            q0 = float(max(theta0[j], 1e-15))
+            hj = min(hj, 0.5 * q0)
+            hj = max(hj, abs_min)
+            th_plus[j] = q0 + hj
+            th_minus[j] = max(q0 - hj, 1e-15)
+
+        g[j] = (F_theta(th_plus) - F_theta(th_minus)) / (2.0 * hj)
+
+    return g
+
+
+def P_with_theta6(theta6, x):
+    """
+    Same posterior risk curve as P_with(), but accepts only the 6 raw parameters:
+    (p, a, b, s, k, th)
+    """
+    p, a, b, s, k, th = np.asarray(theta6, float).ravel()[:6]
+    x = np.asarray(x, float)
+    eps = 1e-12
+
+    logit_val = (
+        np.log(np.clip(p, eps, 1 - eps))
+        - np.log(np.clip(1 - p, eps, 1.0))
+        + dE_full(x, a, b, s, k, th)
+    )
+    return logistic(logit_val)
+
+
+def x50_theta6(theta6, lo=1e-6, hi=6.0, hi_max=100.0):
+    """
+    Solve P(AE|x) = 0.5 using theta6 = (p, a, b, s, k, th),
+    with adaptive bracketing.
+    """
+    f = lambda x: P_with_theta6(theta6, x) - 0.5
+
+    fa = f(lo)
+    fb = f(hi)
+
+    while np.isfinite(fa) and np.isfinite(fb) and fa * fb > 0 and hi < hi_max:
+        hi *= 2.0
+        fb = f(hi)
+
+    if (not np.isfinite(fa)) or (not np.isfinite(fb)) or fa * fb > 0:
+        return np.nan
+
+    try:
+        return float(brentq(f, lo, hi))
+    except Exception:
+        return np.nan
+
+
+def compute_x50_elasticity_rawtheta(
+    theta0,
+    x50_fun=x50_theta6,
+    rel_step=1e-6,
+    abs_min=1e-10,
+    pi_eps=1e-12,
+):
+    """
+    x50 elasticity on RAW theta scale:
+
+      E_x50_j = | (theta_j / x50) * d x50 / d theta_j |
+
+    where theta0 = [p, a, b, s, k, th]
+    """
+    theta0 = np.asarray(theta0, float).ravel()
+    if theta0.size != 6:
+        raise ValueError(f"Expected theta0 of length 6, got {theta0.size}")
+
+    x50_0 = float(x50_fun(theta0))
+
+    d_x_dth = grad_central_theta(
+        x50_fun, theta0, rel_step=rel_step, abs_min=abs_min, pi_eps=pi_eps
+    )
+
+    eps = 1e-12
+    x_safe = max(abs(x50_0), eps)
+
+    th_safe = theta0.copy()
+    th_safe[0] = float(np.clip(th_safe[0], pi_eps, 1 - pi_eps))
+    th_safe[1:] = np.maximum(th_safe[1:], 1e-15)
+
+    elas_x = np.abs(d_x_dth) * np.abs(th_safe) / x_safe
+
+    return {
+        "theta0": theta0,
+        "x50": x50_0,
+        "d_x_dtheta": d_x_dth,
+        "elas_x50": elas_x,
+        "param_names": PARAM_NAMES_X50_ELAS,
+    }
+
+
+def run_x50_elasticity_from_fit_results(
+    fit_results,
+    rel_step=1e-6,
+    abs_min=1e-10,
+    pi_eps=1e-12,
+):
+    """
+    Compute x50 elasticity for FULL (orig) and TRIM directly from fit_results.
+    """
+    elas_orig = compute_x50_elasticity_rawtheta(
+        theta6_from_hat(fit_results["theta_orig"]),
+        x50_fun=x50_theta6,
+        rel_step=rel_step,
+        abs_min=abs_min,
+        pi_eps=pi_eps,
+    )
+
+    elas_trim = compute_x50_elasticity_rawtheta(
+        theta6_from_hat(fit_results["theta_trim"]),
+        x50_fun=x50_theta6,
+        rel_step=rel_step,
+        abs_min=abs_min,
+        pi_eps=pi_eps,
+    )
+
+    return {
+        "orig": elas_orig,
+        "trim": elas_trim,
+    }
+
+
+def make_x50_elasticity_table(elas_res, dataset_name="FULL"):
+    """
+    Tidy x50 elasticity table.
+    """
+    rows = []
+    for name, ex in zip(
+        elas_res["param_names"],
+        elas_res["elas_x50"],
+    ):
+        rows.append(
+            {
+                "Dataset": dataset_name,
+                "Parameter": name,
+                "Elasticity_x50": float(ex),
+            }
+        )
+    return pd.DataFrame(rows)
+
+
+def plot_x50_elasticity_bars(elas_full, elas_trim, figsize=(7, 5), dpi=150):
+    """
+    One-panel bar plot for x50 elasticity.
+    """
+    names = elas_full["param_names"]
+    x = np.arange(len(names))
+    width = 0.36
+
+    fig, ax = plt.subplots(figsize=figsize, dpi=dpi)
+
+    ax.bar(
+        x - width / 2,
+        elas_full["elas_x50"],
+        width=width,
+        label="FULL",
+        alpha=0.85,
+    )
+    ax.bar(
+        x + width / 2,
+        elas_trim["elas_x50"],
+        width=width,
+        label="TRIM",
+        alpha=0.85,
+    )
+
+    ax.set_xticks(x)
+    ax.set_xticklabels(names)
+    ax.set_ylabel("Elasticity", fontweight="bold")
+    ax.set_title("Elasticity of x50", fontweight="bold")
+    ax.grid(axis="y", alpha=0.25)
+    ax.legend(frameon=False)
+
+    plt.tight_layout()
+    return fig, ax

bayesian_noise.py

@@ -0,0 +1,402 @@
+# bayesian_noise.py
+# ============================================================
+# Bayesian noise propagation on RAW X
+# Reference location and Δ widths are evaluated at exact clean-fit x50
+# Function-only module for notebook use
+# ============================================================
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import brentq
+
+from bayesian import load_xy, fit_bayes, P_with, make_trimmed_dataset
+
+# ---------------- Defaults ----------------
+XMAX = 5.0
+GRID_N = 1000
+
+SIGMA_MULT = 0.129
+SIGMA_ADD = 0.144
+
+N_REFIT = 150
+N_TTA = 300
+SEED = 1234
+
+CI_LEVEL = 0.95
+ALPHA = 1.0 - CI_LEVEL
+Q_LO, Q_MD, Q_HI = ALPHA / 2.0, 0.5, 1.0 - ALPHA / 2.0
+
+USE_PRIOR_P = True
+PRIOR_R = (1.05, 1.05)
+TAU = 25.0
+
+plt.rcParams["legend.frameon"] = False
+plt.rcParams["axes.titleweight"] = "bold"
+plt.rcParams["axes.labelweight"] = "bold"
+
+
+# ============================================================
+# 1) DATA
+# ============================================================
+
+def load_original(
+    perc=95,
+    suv_path="suv_percentilesSLOthenUWM.mat",
+    flags_path="flags_combined.mat",
+):
+    return load_xy(perc=perc, suv_path=suv_path, flags_path=flags_path)
+
+
+def load_trimmed(
+    perc=95,
+    suv_path="suv_percentilesSLOthenUWM.mat",
+    flags_path="flags_combined.mat",
+    value_to_drop=2.48122597,
+    tol=1e-3,
+):
+    X, y = load_xy(perc=perc, suv_path=suv_path, flags_path=flags_path)
+    ds = make_trimmed_dataset(X, y, value_to_drop=value_to_drop, tol=tol)
+    return ds["X_trim"], ds["y_trim"]
+
+
+# ============================================================
+# 2) FIT + X50
+# ============================================================
+
+def fit_once(X, y, rng, use_prior_p=USE_PRIOR_P, prior_r=PRIOR_R, tau=TAU):
+    theta, res = fit_bayes(
+        X,
+        y,
+        seed=int(rng.integers(0, 10_000_000)),
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+    return theta, res
+
+
+def x_at_p(theta_hat, p_target=0.5, lo=1e-6, hi=6.0, hi_max=100.0):
+    """
+    Exact root solve for P(AE|x)=p_target.
+    """
+    f = lambda x: P_with(theta_hat, x) - p_target
+    fa, fb = f(lo), f(hi)
+
+    while np.isfinite(fa) and np.isfinite(fb) and fa * fb > 0 and hi < hi_max:
+        hi *= 2.0
+        fb = f(hi)
+
+    if (not np.isfinite(fa)) or (not np.isfinite(fb)) or fa * fb > 0:
+        return np.nan
+
+    try:
+        return float(brentq(f, lo, hi))
+    except Exception:
+        return np.nan
+
+
+# ============================================================
+# 3) NOISE
+# ============================================================
+
+def add_noise_mult(X, sigma, rng):
+    X = np.asarray(X, float)
+    return np.clip(X * np.exp(rng.normal(0, sigma, size=X.shape)), 1e-12, None)
+
+
+def add_noise_add(X, sigma, rng):
+    X = np.asarray(X, float)
+    return np.clip(X + rng.normal(0, sigma, size=X.shape), 1e-12, None)
+
+
+# ============================================================
+# 4) BAND HELPERS
+# ============================================================
+
+def monotone_nondec(y):
+    y = np.asarray(y, float)
+    return np.maximum.accumulate(np.clip(y, 0, 1))
+
+
+def band_quantiles(curves, q_lo=Q_LO, q_md=Q_MD, q_hi=Q_HI):
+    C = np.vstack(curves)
+    ql = monotone_nondec(np.quantile(C, q_lo, axis=0))
+    qm = monotone_nondec(np.quantile(C, q_md, axis=0))
+    qh = monotone_nondec(np.quantile(C, q_hi, axis=0))
+    return ql, qm, qh
+
+
+def delta_width_at_x(lo, hi, xc, x0):
+    if not np.isfinite(x0):
+        return np.nan
+    lo_x = float(np.interp(x0, xc, lo))
+    hi_x = float(np.interp(x0, xc, hi))
+    return hi_x - lo_x
+
+
+# ============================================================
+# 5) CORE ANALYSIS
+# ============================================================
+
+def build_bands_for_dataset(
+    X,
+    y,
+    sigma_mult=SIGMA_MULT,
+    sigma_add=SIGMA_ADD,
+    x_max=XMAX,
+    grid_n=GRID_N,
+    n_refit=N_REFIT,
+    n_tta=N_TTA,
+    seed=SEED,
+    use_prior_p=USE_PRIOR_P,
+    prior_r=PRIOR_R,
+    tau=TAU,
+):
+    rng = np.random.default_rng(seed)
+    xc = np.linspace(1e-12, x_max, grid_n)
+
+    theta_clean, res_clean = fit_once(
+        X, y, rng,
+        use_prior_p=use_prior_p,
+        prior_r=prior_r,
+        tau=tau,
+    )
+    clean_curve = P_with(theta_clean, xc)
+    x50 = x_at_p(theta_clean, p_target=0.5, lo=1e-6, hi=max(6.0, x_max))
+
+    # multiplicative: refit
+    curves_refit_m = []
+    for _ in range(n_refit):
+        Xn = add_noise_mult(X, sigma_mult, rng)
+        thetab, resb = fit_once(Xn, y, rng, use_prior_p=use_prior_p, prior_r=prior_r, tau=tau)
+        if resb.success and np.isfinite(resb.fun):
+            curves_refit_m.append(P_with(thetab, xc))
+    if not curves_refit_m:
+        raise RuntimeError("No successful multiplicative-refit curves.")
+    lo_rm, md_rm, hi_rm = band_quantiles(curves_refit_m)
+
+    # multiplicative: TTA
+    curves_tta_m = []
+    for _ in range(n_tta):
+        xc_n = np.clip(xc * np.exp(rng.normal(0, sigma_mult, size=xc.size)), 1e-12, None)
+        curves_tta_m.append(P_with(theta_clean, xc_n))
+    lo_tm, md_tm, hi_tm = band_quantiles(curves_tta_m)
+
+    # additive: refit
+    curves_refit_a = []
+    for _ in range(n_refit):
+        Xn = add_noise_add(X, sigma_add, rng)
+        thetab, resb = fit_once(Xn, y, rng, use_prior_p=use_prior_p, prior_r=prior_r, tau=tau)
+        if resb.success and np.isfinite(resb.fun):
+            curves_refit_a.append(P_with(thetab, xc))
+    if not curves_refit_a:
+        raise RuntimeError("No successful additive-refit curves.")
+    lo_ra, md_ra, hi_ra = band_quantiles(curves_refit_a)
+
+    # additive: TTA
+    curves_tta_a = []
+    for _ in range(n_tta):
+        xc_n = np.clip(xc + rng.normal(0, sigma_add, size=xc.size), 1e-12, None)
+        curves_tta_a.append(P_with(theta_clean, xc_n))
+    lo_ta, md_ta, hi_ta = band_quantiles(curves_tta_a)
+
+    return {
+        "xc": xc,
+        "theta_clean": theta_clean,
+        "fit_success": bool(res_clean.success),
+        "clean_curve": clean_curve,
+        "x50": x50,
+        "mult": {
+            "sigma": sigma_mult,
+            "refit": (lo_rm, md_rm, hi_rm),
+            "tta": (lo_tm, md_tm, hi_tm),
+            "D_refit": delta_width_at_x(lo_rm, hi_rm, xc, x50),
+            "D_tta": delta_width_at_x(lo_tm, hi_tm, xc, x50),
+            "n_refit_success": len(curves_refit_m),
+            "n_tta": len(curves_tta_m),
+        },
+        "add": {
+            "sigma": sigma_add,
+            "refit": (lo_ra, md_ra, hi_ra),
+            "tta": (lo_ta, md_ta, hi_ta),
+            "D_refit": delta_width_at_x(lo_ra, hi_ra, xc, x50),
+            "D_tta": delta_width_at_x(lo_ta, hi_ta, xc, x50),
+            "n_refit_success": len(curves_refit_a),
+            "n_tta": len(curves_tta_a),
+        },
+    }
+
+
+def run_noise_analysis(
+    perc=95,
+    suv_path="suv_percentilesSLOthenUWM.mat",
+    flags_path="flags_combined.mat",
+    value_to_drop=2.48122597,
+    tol=1e-3,
+    sigma_mult=SIGMA_MULT,
+    sigma_add=SIGMA_ADD,
+    x_max=XMAX,
+    grid_n=GRID_N,
+    n_refit=N_REFIT,
+    n_tta=N_TTA,
+    seed=SEED,
+    use_prior_p=USE_PRIOR_P,
+    prior_r=PRIOR_R,
+    tau=TAU,
+):
+    X_full, y_full = load_original(
+        perc=perc, suv_path=suv_path, flags_path=flags_path
+    )
+    X_trim, y_trim = load_trimmed(
+        perc=perc, suv_path=suv_path, flags_path=flags_path,
+        value_to_drop=value_to_drop, tol=tol
+    )
+
+    full_res = build_bands_for_dataset(
+        X_full, y_full,
+        sigma_mult=sigma_mult, sigma_add=sigma_add,
+        x_max=x_max, grid_n=grid_n,
+        n_refit=n_refit, n_tta=n_tta,
+        seed=seed,
+        use_prior_p=use_prior_p, prior_r=prior_r, tau=tau,
+    )
+    trim_res = build_bands_for_dataset(
+        X_trim, y_trim,
+        sigma_mult=sigma_mult, sigma_add=sigma_add,
+        x_max=x_max, grid_n=grid_n,
+        n_refit=n_refit, n_tta=n_tta,
+        seed=seed + 100,
+        use_prior_p=use_prior_p, prior_r=prior_r, tau=tau,
+    )
+
+    return {"full": full_res, "trim": trim_res}
+
+
+# ============================================================
+# 6) SUMMARY TABLE
+# ============================================================
+
+def make_noise_summary_table(noise_res):
+    rows = []
+    for dataset_name, block in [("FULL", noise_res["full"]), ("TRIM", noise_res["trim"])]:
+        for noise_name, key in [("multiplicative", "mult"), ("additive", "add")]:
+            rows.append({
+                "Dataset": dataset_name,
+                "NoiseType": noise_name,
+                "Sigma": float(block[key]["sigma"]),
+                "x50": float(block["x50"]),
+                "Delta_refit_at_x50": float(block[key]["D_refit"]),
+                "Delta_tta_at_x50": float(block[key]["D_tta"]),
+                "N_refit_success": int(block[key]["n_refit_success"]),
+                "N_tta": int(block[key]["n_tta"]),
+            })
+    import pandas as pd
+    return pd.DataFrame(rows)
+
+
+# ============================================================
+# 7) PLOTTING
+# ============================================================
+
+def plot_4panels(
+    full_res,
+    trim_res,
+    save_prefix="noise_cb_raw_4panels_sigma129_0144_x50",
+    x_max=XMAX,
+    sigma_mult=SIGMA_MULT,
+    sigma_add=SIGMA_ADD,
+):
+    fig, axs = plt.subplots(2, 2, figsize=(12.5, 7.5), dpi=160, sharex=True, sharey=True)
+
+    COL = {"mult": "#1f78b4", "add": "#ff7f00"}
+    ALP = {"refit": 0.25, "tta": 0.12}
+
+    def draw_panel(ax, res, mode, letter, dataset_label):
+        xc = res["xc"]
+        clean = res["clean_curve"]
+        x50 = float(res["x50"])
+
+        lo_r, _, hi_r = res[mode]["refit"]
+        lo_t, _, hi_t = res[mode]["tta"]
+        Dref = res[mode]["D_refit"]
+        Dtta = res[mode]["D_tta"]
+
+        ax.plot(xc, clean, color="k", lw=2.6, zorder=3)
+        ax.fill_between(xc, lo_r, hi_r, color=COL[mode], alpha=ALP["refit"], zorder=1)
+        ax.fill_between(xc, lo_t, hi_t, color=COL[mode], alpha=ALP["tta"], zorder=0)
+
+        # exact x50 reference
+        if np.isfinite(x50):
+            ax.axvline(x50, color="#666", ls="--", lw=1.2, alpha=0.9, zorder=4)
+
+            # vertical band-width markers at x50
+            lo_r_x = float(np.interp(x50, xc, lo_r))
+            hi_r_x = float(np.interp(x50, xc, hi_r))
+            lo_t_x = float(np.interp(x50, xc, lo_t))
+            hi_t_x = float(np.interp(x50, xc, hi_t))
+
+            ax.plot([x50, x50], [lo_r_x, hi_r_x], color=COL[mode], lw=2.2, alpha=0.95, zorder=6)
+            ax.scatter([x50, x50], [lo_r_x, hi_r_x], color=COL[mode], s=18, alpha=0.95, zorder=7)
+
+            ax.plot([x50, x50], [lo_t_x, hi_t_x], color=COL[mode], lw=2.2, alpha=0.45, zorder=5)
+            ax.scatter([x50, x50], [lo_t_x, hi_t_x], color=COL[mode], s=18, alpha=0.45, zorder=6)
+
+        ax.set_title(letter, fontsize=14, fontweight="bold", pad=8)
+
+        ax.text(
+            0.02, 0.96, dataset_label,
+            transform=ax.transAxes, ha="left", va="top",
+            fontsize=10, color="#111"
+        )
+        ax.text(
+            0.02, 0.88, f"x50={x50:.2f}",
+            transform=ax.transAxes, ha="left", va="top",
+            fontsize=10, color="#111"
+        )
+        ax.text(
+            0.03, 0.10, f"Δr@x50={Dref:.2f}  Δt@x50={Dtta:.2f}",
+            transform=ax.transAxes, ha="left", va="center", fontsize=9,
+            bbox=dict(
+                facecolor="white",
+                edgecolor=COL[mode],
+                boxstyle="round,pad=0.25,rounding_size=0.02",
+                lw=0.9, alpha=0.95
+            )
+        )
+
+        ax.set_xlim(0, x_max)
+        ax.set_ylim(-0.05, 1.05)
+        ax.grid(alpha=0.25)
+
+    draw_panel(axs[0, 0], full_res, "mult", "A", "FULL (F)")
+    draw_panel(axs[0, 1], full_res, "add",  "B", "FULL (F)")
+    draw_panel(axs[1, 0], trim_res, "mult", "C", "TRIM (T)")
+    draw_panel(axs[1, 1], trim_res, "add",  "D", "TRIM (T)")
+
+    axs[1, 0].set_xlabel("X")
+    axs[1, 1].set_xlabel("X")
+    axs[0, 0].set_ylabel("P(AE | x)")
+    axs[1, 0].set_ylabel("P(AE | x)")
+
+    handles = [
+        plt.Line2D([0], [0], color="k", lw=2.6, label="clean Bayesian fit"),
+        plt.Rectangle((0, 0), 1, 1, facecolor=COL["mult"], alpha=ALP["refit"],
+                      label=f"refit band, mult σ={sigma_mult:.3f}"),
+        plt.Rectangle((0, 0), 1, 1, facecolor=COL["mult"], alpha=ALP["tta"],
+                      label=f"TTA band, mult σ={sigma_mult:.3f}"),
+        plt.Rectangle((0, 0), 1, 1, facecolor=COL["add"], alpha=ALP["refit"],
+                      label=f"refit band, add σ={sigma_add:.3f}"),
+        plt.Rectangle((0, 0), 1, 1, facecolor=COL["add"], alpha=ALP["tta"],
+                      label=f"TTA band, add σ={sigma_add:.3f}"),
+        plt.Line2D([0], [0], color="#666", lw=1.2, ls="--", label="x50"),
+    ]
+    fig.legend(handles, [h.get_label() for h in handles],
+               loc="lower center", ncol=3, fontsize=10)
+
+    plt.tight_layout(rect=[0, 0.08, 1, 1])
+
+    if save_prefix is not None:
+        fig.savefig(f"{save_prefix}.png", dpi=300, bbox_inches="tight")
+        fig.savefig(f"{save_prefix}.pdf", bbox_inches="tight")
+
+    return fig, axs

bayesianconstraints.py

@@ -0,0 +1,139 @@
+# bayesianconstraints.py
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+
+def logistic(z):
+    z = np.clip(z, -60, 60)
+    return 1.0 / (1.0 + np.exp(-z))
+
+
+def logit(p):
+    p = np.clip(p, 1e-12, 1 - 1e-12)
+    return np.log(p) - np.log(1 - p)
+
+
+def dE_gamma(x, a, b, k, theta, s=1.0):
+    x = np.asarray(x, float)
+    return (a - k) * np.log(x) - (a + b) * np.log1p(x / s) + x / theta
+
+
+def dE_lognorm(x, a, b, mu, sigma, s=1.0):
+    x = np.asarray(x, float)
+    return ((mu - np.log(x)) ** 2) / (2.0 * sigma ** 2) + a * np.log(x) - (a + b) * np.log1p(x / s)
+
+
+def check_gamma_constraints(a, b, k, theta, s):
+    positivity = (a > 0) and (b > 0) and (k > 0) and (theta > 0) and (s > 0)
+
+    if not positivity:
+        status = "FAIL"
+        note = "all parameters must be > 0"
+    elif a > k:
+        status = "MEET"
+        note = "a > k"
+    elif a == k:
+        status = "EDGE"
+        note = "a = k"
+    else:
+        status = "FAIL"
+        note = "a < k"
+
+    return {
+        "model": "BetaPrime vs Gamma",
+        "a": a,
+        "b": b,
+        "k": k,
+        "theta": theta,
+        "s": s,
+        "all_positive": positivity,
+        "a_ge_k": positivity and (a >= k),
+        "status": status,
+        "note": note,
+    }
+
+
+def check_lognorm_constraints(a, b, mu, sigma, s):
+    positivity = (a > 0) and (b > 0) and (sigma > 0) and (s > 0)
+
+    if not positivity:
+        status = "FAIL"
+        note = "a,b,sigma,s must be > 0"
+    else:
+        status = "FAIL"
+        note = "left tail goes to 1, not 0"
+
+    return {
+        "model": "BetaPrime vs LogNormal",
+        "a": a,
+        "b": b,
+        "mu": mu,
+        "sigma": sigma,
+        "s": s,
+        "all_positive": positivity,
+        "status": status,
+        "note": note,
+    }
+
+
+def make_gamma_constraint_table(param_list):
+    return pd.DataFrame([check_gamma_constraints(**p) for p in param_list])
+
+
+def make_lognorm_constraint_table(param_list):
+    return pd.DataFrame([check_lognorm_constraints(**p) for p in param_list])
+
+
+def plot_gamma_constraint_curves(param_meet, param_edge_fail, p_prior=5/58,
+                                 x=None, figsize=(10, 6)):
+    if x is None:
+        x = np.logspace(-12, 3, 900)
+
+    c0 = logit(p_prior)
+    fig, ax = plt.subplots(figsize=figsize)
+
+    for P in param_meet:
+        p = logistic(c0 + dE_gamma(x, **P))
+        label = f"MEET (a={P['a']}, b={P['b']}, k={P['k']}, θ={P['theta']}, s={P['s']})"
+        ax.plot(x, p, label=label)
+
+    for P in param_edge_fail:
+        p = logistic(c0 + dE_gamma(x, **P))
+        tag = "EDGE" if P["a"] == P["k"] else "FAIL"
+        label = f"{tag} (a={P['a']}, b={P['b']}, k={P['k']}, θ={P['theta']}, s={P['s']})"
+        ax.plot(x, p, linestyle="--", label=label)
+
+    ax.set_xscale("log")
+    ax.set_ylim(-0.05, 1.05)
+    ax.set_xlabel("x")
+    ax.set_ylabel("P(AE | x)")
+    ax.set_title("Gamma NC — posterior-like curves")
+    ax.grid(True, which="both")
+    ax.legend(fontsize=8)
+    fig.tight_layout()
+    return fig, ax
+
+
+def plot_lognorm_constraint_curves(param_list, p_prior=5/58,
+                                   x=None, figsize=(10, 6)):
+    if x is None:
+        x = np.logspace(-12, 3, 900)
+
+    c0 = logit(p_prior)
+    fig, ax = plt.subplots(figsize=figsize)
+
+    for P in param_list:
+        p = logistic(c0 + dE_lognorm(x, **P))
+        label = f"(a={P['a']}, b={P['b']}, μ={P['mu']}, σ={P['sigma']}, s={P['s']})"
+        ax.plot(x, p, label=label)
+
+    ax.set_xscale("log")
+    ax.set_ylim(-0.05, 1.05)
+    ax.set_xlabel("x")
+    ax.set_ylabel("P(AE | x)")
+    ax.set_title("LogNormal NC — posterior-like curves")
+    ax.grid(True, which="both")
+    ax.legend(fontsize=8)
+    fig.tight_layout()
+    return fig, ax