Delete 'src/irae_risk/full_logistic_Zahra'

src/irae_risk/full_logistic_Zahra

@@ -1,2326 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from scipy.optimize import minimize
-
-
-# ============================================================
-# Stable sigmoid
-# ============================================================
-
-def _sigmoid_stable(z):
-    z = np.asarray(z, float)
-    z = np.clip(z, -50.0, 50.0)
-    return 1.0 / (1.0 + np.exp(-z))
-
-
-# ============================================================
-# 1) Model
-# ============================================================
-
-def model_p(x, b):
-    """p(x|b) = sigmoid(b0 + b1*x)."""
-    x = np.asarray(x, float).reshape(-1)
-    b0, b1 = np.asarray(b, float).reshape(2)
-    return _sigmoid_stable(b0 + b1 * x)
-
-
-def design_matrix(x):
-    """Design matrix X = [1, x]."""
-    x = np.asarray(x, float).reshape(-1)
-    return np.column_stack([np.ones_like(x), x])
-
-
-# ============================================================
-# 2) Likelihood
-# ============================================================
-
-def nll(x, y, b, l2=0.0):
-    """
-    Penalized negative log-likelihood:
-        NLL(b) = -sum[y log p + (1-y) log(1-p)] + 0.5*l2*||b||^2
-    """
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, float).reshape(-1)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x, b)
-    eps = 1e-12
-    p = np.clip(p, eps, 1 - eps)
-
-    base = -np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))
-    pen = 0.5 * l2 * float(np.dot(b, b))
-    return base + pen
-
-
-def llf(x, y, b):
-    """
-    Ordinary (unpenalized) log-likelihood at fitted parameters.
-    """
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, float).reshape(-1)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x, b)
-    eps = 1e-12
-    p = np.clip(p, eps, 1 - eps)
-    return float(np.sum(y * np.log(p) + (1 - y) * np.log(1 - p)))
-
-
-# ============================================================
-# 3) Gradient / Hessian / Covariance
-# ============================================================
-
-def grad_nll(x, y, b, l2=0.0):
-    """
-    Gradient of penalized NLL:
-        g(b) = X^T (p - y) + l2*b
-    """
-    X = design_matrix(x)
-    y = np.asarray(y, float).reshape(-1)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x, b)
-    return X.T @ (p - y) + l2 * b
-
-
-def hess_nll(x, b, l2=0.0):
-    """
-    Hessian of penalized NLL:
-        H(b) = X^T W X + l2*I
-        W = diag(p*(1-p))
-    """
-    X = design_matrix(x)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x, b)
-    w = p * (1 - p)
-    return X.T @ (w[:, None] * X) + l2 * np.eye(2)
-
-
-def covariance(x, b, l2=0.0):
-    """
-    Cov(b) ≈ H(b)^(-1), where H is the penalized Hessian if l2 > 0.
-    Robust to near-singular Hessians.
-    """
-    H = hess_nll(x, b, l2=l2)
-    try:
-        return np.linalg.inv(H)
-    except np.linalg.LinAlgError:
-        return np.linalg.pinv(H)
-
-
-def standard_errors(x, b, l2=0.0):
-    """
-    SE = sqrt(diag(Cov)).
-    """
-    C = covariance(x, b, l2=l2)
-    return np.sqrt(np.maximum(np.diag(C), 0.0))
-
-
-# Compatibility alias
-def logit_poly_cov(x, b, l2=0.0):
-    return covariance(x, b, l2=l2)
-
-
-# ============================================================
-# 4) Fit
-# ============================================================
-def fit_newton(x, y, b_start=None, max_iter=50, tol=1e-8, l2=0.0):
-    """
-    Newton updates for penalized NLL with backtracking line-search.
-
-    Update:
-        b_new = b - alpha * H^{-1} g
-    alpha shrinks until NLL decreases.
-    """
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-
-    if b_start is None:
-        b = np.array([0.0, 0.0], float)
-    else:
-        b = np.asarray(b_start, float).reshape(2)
-
-    f = nll(x, y, b, l2=l2)
-
-    for _ in range(max_iter):
-        g = grad_nll(x, y, b, l2=l2)
-        H = hess_nll(x, b, l2=l2)
-
-        try:
-            step = np.linalg.solve(H, g)
-        except np.linalg.LinAlgError:
-            step = np.linalg.pinv(H) @ g
-
-        alpha = 1.0
-        while alpha > 1e-6:
-            b_new = b - alpha * step
-            f_new = nll(x, y, b_new, l2=l2)
-            if np.isfinite(f_new) and f_new <= f:
-                break
-            alpha *= 0.5
-
-        if alpha <= 1e-6:
-            break
-
-        if np.max(np.abs(b_new - b)) < tol:
-            b = b_new
-            break
-
-        b, f = b_new, f_new
-
-    return b
-
-
-# ============================================================
-# 14) Overlay plot (LOG left, RAW right)
-# ============================================================
-def plot_overlay_two_panels_final(
-    r_log_full, r_log_trim, r_raw_full, r_raw_trim,
-    dy_full=-0.010, dy_trim=0.010
-):
-    import numpy as np
-    import matplotlib.pyplot as plt
-
-    COL_NC   = "#4C78A8"
-    COL_AE   = "#F58518"
-    COL_FULL = "#1f77b4"
-    COL_TRIM = "#ff7f0e"
-
-    fig, axes = plt.subplots(1, 2, figsize=(9.5, 3.5), sharey=True)
-    ax1, ax2 = axes
-
-    def draw_panel(ax, r_full, r_trim, xlabel, panel_label,
-                   show_legend=False):
-
-        xF = np.asarray(r_full["x"], float)
-        yF = np.asarray(r_full["y"], int)
-
-        xT = np.asarray(r_trim["x"], float)
-        yT = np.asarray(r_trim["y"], int)
-
-        xx = np.linspace(
-            min(xF.min(), xT.min()),
-            max(xF.max(), xT.max()),
-            500
-        )
-
-        # keep x-values unchanged
-        xF_plot = xF
-        xT_plot = xT
-
-        # vertical offsets only
-        yF_plot = yF + np.where(yF == 0, dy_full, -dy_full)
-        yT_plot = yT + np.where(yT == 0, dy_trim, -dy_trim)
-
-        # FULL = filled markers
-        ax.scatter(
-            xF_plot[yF == 0], yF_plot[yF == 0],
-            s=16,
-            color=COL_NC,
-            alpha=0.70,
-            edgecolors="none",
-            label="data: FULL NC",
-            zorder=3
-        )
-
-        ax.scatter(
-            xF_plot[yF == 1], yF_plot[yF == 1],
-            s=16,
-            color=COL_AE,
-            alpha=0.80,
-            edgecolors="none",
-            label="data: FULL AE",
-            zorder=3
-        )
-
-        # TRIM = outlined markers
-        ax.scatter(
-            xT_plot[yT == 0], yT_plot[yT == 0],
-            s=24,
-            facecolors=COL_NC,
-            edgecolors="black",
-            linewidths=0.45,
-            alpha=0.95,
-            label="data: TRIM NC",
-            zorder=4
-        )
-
-        ax.scatter(
-            xT_plot[yT == 1], yT_plot[yT == 1],
-            s=24,
-            facecolors=COL_AE,
-            edgecolors="black",
-            linewidths=0.45,
-            alpha=0.95,
-            label="data: TRIM AE",
-            zorder=4
-        )
-
-        # logistic fits
-        ax.plot(
-            xx,
-            model_p(xx, r_full["b"]),
-            lw=1.8,
-            color=COL_FULL,
-            label="fit FULL",
-            zorder=2
-        )
-
-        ax.plot(
-            xx,
-            model_p(xx, r_trim["b"]),
-            lw=1.8,
-            color=COL_TRIM,
-            label="fit TRIM",
-            zorder=2
-        )
-
-        # legend only in panel B
-        if show_legend:
-            ax.legend(
-                loc="lower right",
-                fontsize=7,
-                markerscale=0.9,
-                frameon=True,
-                framealpha=1.0,
-                edgecolor="0.7",
-                handlelength=1.8,
-                borderpad=0.4,
-                labelspacing=0.4,
-                handletextpad=0.5
-            )
-
-        ax.set_xlabel(xlabel)
-        ax.set_ylim(-0.08, 1.08)
-
-        ax.text(
-            0.05, 0.90,
-            panel_label,
-            transform=ax.transAxes,
-            fontsize=11
-        )
-
-    draw_panel(
-        ax1,
-        r_log_full,
-        r_log_trim,
-        "log(X)",
-        "A",
-        show_legend=False
-    )
-
-    draw_panel(
-        ax2,
-        r_raw_full,
-        r_raw_trim,
-        "X",
-        "B",
-        show_legend=True
-    )
-
-    ax1.set_ylabel("P(AE | X = x)")
-
-    for ax in axes:
-        ax.grid(False)
-        ax.tick_params(labelsize=8)
-
-    plt.tight_layout()
-    plt.show()
-
-    return fig, axes
-
-    return fig, axes
-# ============================================================
-# 5) Goodness of fit
-# ============================================================
-
-def goodness_of_fit(x, y, b, thresh=0.5, l2=0.0):
-    """
-    Returns:
-      LLF, NLL, AIC, BIC, Accuracy, n, k
-
-    Notes
-    -----
-    Fit may use l2 > 0, but GOF metrics below are computed from the
-    ordinary (unpenalized) likelihood evaluated at the fitted parameters.
-    The argument l2 is kept only for interface consistency.
-    """
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x, b)
-    eps = 1e-12
-    p = np.clip(p, eps, 1 - eps)
-
-    LLF = np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))
-    NLL = -LLF
-
-    n = len(x)
-    k = len(b)
-
-    AIC = 2 * k - 2 * LLF
-    BIC = k * np.log(n) - 2 * LLF
-
-    yhat = (p >= thresh).astype(int)
-    acc = np.mean(yhat == y)
-
-    return {
-        "LLF": float(LLF),
-        "NLL": float(NLL),
-        "AIC": float(AIC),
-        "BIC": float(BIC),
-        "A": float(acc),
-        "n": int(n),
-        "k": int(k),
-    }
-
-
-# ============================================================
-# 6) x50 / Wald helpers / compact fit
-# ============================================================
-def x50(b):
-    """
-    Model-scale midpoint:
-        x50 = -b0 / b1
-
-    For LOG panels, this is on the log(x) scale.
-    Raw-scale SUV50 is exp(x50).
-    """
-    b0, b1 = np.asarray(b, float).reshape(2)
-    return np.nan if np.abs(b1) < 1e-12 else (-b0 / b1)
-
-
-def check_x50_consistency(P):
-    """
-    Diagnostic check for x50 consistency.
-
-    For LOG models:
-        x50_model is on log(X) scale
-        SUV50 is on raw X scale = exp(x50_model)
-
-    For RAW models:
-        x50_model = SUV50
-
-    Correct result:
-        P(x50_model) should be approximately 0.5
-    """
-    for key, pk in P.items():
-        b = np.asarray(pk["b"], float).reshape(2)
-        trans = pk.get("transform", "")
-
-        x50_model = x50(b)
-        suv50 = np.exp(x50_model) if trans == "log" else x50_model
-        p_at_x50 = model_p(np.array([x50_model]), b)[0]
-
-        print(
-            key,
-            "| transform =", trans,
-            "| x50_model =", x50_model,
-            "| SUV50 =", suv50,
-            "| P(x50) =", p_at_x50
-        )
-
-
-def x50_wald_ci(b, cov, z=1.959963984540054):
-    """
-    Wald CI for x50 = -b0/b1 via delta method.
-    Returned on MODEL scale.
-    """
-    b = np.asarray(b, float).reshape(2)
-    cov = np.asarray(cov, float).reshape(2, 2)
-    b0, b1 = b
-
-    if np.abs(b1) < 1e-12:
-        return np.nan, np.nan
-
-    xhat = -b0 / b1
-    grad = np.array([-1.0 / b1, b0 / (b1 ** 2)], float)
-    var = float(grad.T @ cov @ grad)
-    se = np.sqrt(max(var, 0.0))
-    return float(xhat - z * se), float(xhat + z * se)
-
-
-def wald_ci(b, cov, z=1.959963984540054):
-    """
-    Wald CI for parameters: b_i ± z*SE_i.
-    """
-    b = np.asarray(b, float).reshape(2)
-    cov = np.asarray(cov, float).reshape(2, 2)
-    se = np.sqrt(np.maximum(np.diag(cov), 0.0))
-    return b - z * se, b + z * se
-
-
-def fit_pack(x, y, name="", thresh=0.5, l2=0.0, z=1.959963984540054):
-    """
-    Fit + covariance + GOF + parameter Wald CI.
-    x should already be on the MODEL scale.
-    """
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-
-    b = fit_newton(x, y, l2=l2)
-    cov = covariance(x, b, l2=l2)
-    gof = goodness_of_fit(x, y, b, thresh=thresh, l2=l2)
-    lcl, ucl = wald_ci(b, cov, z=z)
-
-    return {
-        "name": name,
-        "x": x,
-        "y": y,
-        "b": b,
-        "cov": cov,
-        "gof": gof,
-        "LCL": lcl,
-        "UCL": ucl,
-        "l2": float(l2),
-    }
-
-
-def trim_nc_by_value(x_raw, y, target=2.48, tol=0.05):
-    """
-    Remove ONE NC sample (y==0) with x_raw closest to target.
-    """
-    x_raw = np.asarray(x_raw, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-
-    nc_idx = np.where(y == 0)[0]
-    if len(nc_idx) == 0:
-        raise ValueError("No NC samples found (y==0).")
-
-    j = nc_idx[np.argmin(np.abs(x_raw[nc_idx] - target))]
-    diff = float(np.abs(x_raw[j] - target))
-    if diff > tol:
-        print(f"[trim warning] closest NC to {target} is {x_raw[j]:.6f} (diff={diff:.6f}) > tol={tol}")
-
-    mask = np.ones_like(y, dtype=bool)
-    mask[j] = False
-    print(f"[trim] removed index={j}, x_raw={x_raw[j]:.6f}, y={y[j]}")
-    return x_raw[mask], y[mask]
-
-def s50_from_b(b, transform="raw"):
-    """
-    Local slope s50 = dp/dx at x50 on RAW x scale.
-    """
-    b = np.asarray(b, float).reshape(2)
-    b0, b1 = b
-
-    if np.abs(b1) < 1e-12:
-        return np.nan
-
-    x50_model = x50(b)
-
-    if transform == "raw":
-        return float(b1 / 4.0)
-
-    elif transform == "log":
-        x50_raw = np.exp(x50_model)
-        return float(b1 / (4.0 * x50_raw))
-
-    else:
-        raise ValueError("transform must be 'raw' or 'log'")
-
-
-def s50_normal_ci_from_mvnorm(
-    b, cov, transform="raw",
-    M=200000, seed=123, alpha=0.05,
-    enforce_positive_slope=True, slope_eps=1e-10
-):
-    """
-    Normal-on-MLE CI for s50.
-    """
-    rng = np.random.default_rng(seed)
-
-    b = np.asarray(b, float).reshape(2)
-    cov = np.asarray(cov, float).reshape(2, 2)
-
-    vals = []
-    tries = 0
-    max_tries = 20 * M
-
-    while len(vals) < M and tries < max_tries:
-        tries += 1
-
-        bb = rng.multivariate_normal(mean=b, cov=cov)
-
-        if not np.all(np.isfinite(bb)):
-            continue
-
-        if enforce_positive_slope and bb[1] <= slope_eps:
-            continue
-
-        val = s50_from_b(bb, transform=transform)
-
-        if np.isfinite(val):
-            vals.append(val)
-
-    if len(vals) == 0:
-        return np.nan, np.nan, np.nan, 0
-
-    vals = np.asarray(vals, float)
-    q = np.quantile(vals, [alpha / 2, 0.5, 1.0 - alpha / 2])
-
-    return float(q[0]), float(q[1]), float(q[2]), int(len(vals))
-
-
-def s50_wald_ci_numeric(
-    b, cov, transform="raw",
-    z=1.959963984540054,
-    eps=1e-5
-):
-    """
-    Delta-method CI for s50 using numerical derivatives.
-    """
-    b = np.asarray(b, float).reshape(2)
-    cov = np.asarray(cov, float).reshape(2, 2)
-
-    s_hat = s50_from_b(b, transform=transform)
-
-    if not np.isfinite(s_hat):
-        return np.nan, np.nan
-
-    grad = np.zeros(2, float)
-
-    for j in range(2):
-        step = eps * max(1.0, abs(b[j]))
-
-        bp = b.copy()
-        bm = b.copy()
-
-        bp[j] += step
-        bm[j] -= step
-
-        sp = s50_from_b(bp, transform=transform)
-        sm = s50_from_b(bm, transform=transform)
-
-        grad[j] = (sp - sm) / (2.0 * step)
-
-    var = float(grad.T @ cov @ grad)
-    se = np.sqrt(max(var, 0.0))
-
-    return float(s_hat - z * se), float(s_hat + z * se)
-
-
-# ============================================================
-# 7) Alternative confidence-interval estimation
-# ============================================================
-# Final terminology:
-#   Wald          : analytical approximation using the fitted covariance;
-#   MC            : Monte Carlo propagation from the local Gaussian approximation;
-#   Nonparametric : ordinary patient-level nonparametric bootstrap;
-#   Stratified    : class-stratified nonparametric bootstrap, retained for comparison;
-#   Parametric    : model-based Bernoulli bootstrap.
-#
-# The delta method is used internally for analytical propagation under
-# the Wald approximation; it is not treated as a separate method.
-
-from collections import OrderedDict
-
-
-CI_METHODS = [
-    "Wald",
-    "MC",
-    "Nonparametric",
-    "Stratified",
-    "Parametric",
-]
-
-MC_DRAWS_BANDS = 20_000
-MC_DRAWS_TABLE = 200_000
-
-
-def eta_se_grid(x_grid, cov):
-    """Standard error of eta(x) = b0 + b1*x on a model-scale grid."""
-    x_grid = np.asarray(x_grid, float).reshape(-1)
-    cov = np.asarray(cov, float).reshape(2, 2)
-
-    Xg = design_matrix(x_grid)
-    var_eta = np.einsum("ij,jk,ik->i", Xg, cov, Xg)
-    return np.sqrt(np.maximum(var_eta, 0.0))
-
-
-def ci_band_wald(x_grid, b, cov, z=1.959963984540054):
-    """
-    Pointwise Wald confidence band for p(x).
-
-    The fitted-parameter covariance is propagated to the probability scale
-    using the first-order delta method.
-    """
-    x_grid = np.asarray(x_grid, float).reshape(-1)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x_grid, b)
-    se_eta = eta_se_grid(x_grid, cov)
-    se_p = p * (1.0 - p) * se_eta
-
-    lo = np.clip(p - z * se_p, 0.0, 1.0)
-    hi = np.clip(p + z * se_p, 0.0, 1.0)
-    return lo, p, hi
-
-
-# Historical alias retained for notebook compatibility.
-def ci_band_delta(x_grid, b, cov, z=1.959963984540054):
-    return ci_band_wald(x_grid, b, cov, z=z)
-
-
-def gaussian_parameter_draws(
-    b,
-    cov,
-    M=MC_DRAWS_BANDS,
-    seed=123,
-    enforce_positive_slope=True,
-    slope_eps=1e-10,
-    x50_bounds=None,
-):
-    """
-    Draw beta* ~ N(beta_hat, Cov_hat) for Monte Carlo propagation.
-
-    Parameters
-    ----------
-    x50_bounds : tuple(float, float) or None
-        Optional admissible interval for model-scale x50. When supplied,
-        draws with x50 outside [lower, upper] are rejected.
-    """
-    rng = np.random.default_rng(seed)
-
-    b = np.asarray(b, float).reshape(2)
-    cov = np.asarray(cov, float).reshape(2, 2)
-
-    draws = []
-    attempts = 0
-    rejected_nonfinite = 0
-    rejected_slope = 0
-    rejected_x50 = 0
-
-    max_attempts = max(50 * int(M), 1000)
-
-    while len(draws) < int(M) and attempts < max_attempts:
-        attempts += 1
-
-        try:
-            bb = rng.multivariate_normal(mean=b, cov=cov)
-        except Exception:
-            break
-
-        if not np.all(np.isfinite(bb)):
-            rejected_nonfinite += 1
-            continue
-
-        if enforce_positive_slope and bb[1] <= slope_eps:
-            rejected_slope += 1
-            continue
-
-        if x50_bounds is not None:
-            x50_draw = x50(bb)
-
-            if not np.isfinite(x50_draw):
-                rejected_x50 += 1
-                continue
-
-            x50_lower, x50_upper = x50_bounds
-
-            if not (x50_lower <= x50_draw <= x50_upper):
-                rejected_x50 += 1
-                continue
-
-        draws.append(bb)
-
-    arr = np.asarray(draws, float) if draws else np.empty((0, 2), float)
-
-    diagnostics = {
-        "attempted": int(attempts),
-        "successful": int(len(arr)),
-        "rejected": int(attempts - len(arr)),
-        "rejected_nonfinite": int(rejected_nonfinite),
-        "rejected_slope": int(rejected_slope),
-        "rejected_x50": int(rejected_x50),
-        "success_rate": (
-            float(len(arr) / attempts)
-            if attempts > 0 else np.nan
-        ),
-    }
-
-    return arr, diagnostics
-
-def bootstrap_band_from_params(x_grid, pars, alpha=0.05):
-    """Convert parameter draws to pointwise confidence bands."""
-    x_grid = np.asarray(x_grid, float).reshape(-1)
-    pars = np.asarray(pars, float)
-
-    if pars.ndim != 2 or pars.shape[0] == 0:
-        nan = np.full_like(x_grid, np.nan, dtype=float)
-        return nan, nan, nan
-
-    curves = np.asarray([model_p(x_grid, bb) for bb in pars], float)
-    q = np.quantile(curves, [alpha / 2, 0.5, 1.0 - alpha / 2], axis=0)
-    return q[0], q[1], q[2]
-
-
-def ci_band_normal_mle_sim(
-    x_grid,
-    b,
-    cov,
-    M=MC_DRAWS_BANDS,
-    seed=123,
-    alpha=0.05,
-    enforce_positive_slope=True,
-    enforce_x50_in_grid=False,
-    slope_eps=1e-10,
-):
-    """
-    Monte Carlo confidence band from the local Gaussian approximation.
-
-    The historical function name and signature are retained. The old
-    x50-in-grid filter is intentionally ignored.
-    """
-    draws, _ = gaussian_parameter_draws(
-        b,
-        cov,
-        M=M,
-        seed=seed,
-        enforce_positive_slope=enforce_positive_slope,
-        slope_eps=slope_eps,
-    )
-    return bootstrap_band_from_params(x_grid, draws, alpha=alpha)
-
-
-def x50_normal_ci_from_mvnorm(
-    b,
-    cov,
-    M=MC_DRAWS_TABLE,
-    seed=123,
-    alpha=0.05,
-    enforce_positive_slope=True,
-    slope_eps=1e-10,
-):
-    """Monte Carlo interval for model-scale x50 = -b0/b1."""
-    draws, _ = gaussian_parameter_draws(
-        b,
-        cov,
-        M=M,
-        seed=seed,
-        enforce_positive_slope=enforce_positive_slope,
-        slope_eps=slope_eps,
-    )
-
-    if len(draws) == 0:
-        return np.nan, np.nan, np.nan, 0
-
-    vals = np.asarray([x50(bb) for bb in draws], float)
-    vals = vals[np.isfinite(vals)]
-
-    if len(vals) == 0:
-        return np.nan, np.nan, np.nan, 0
-
-    q = np.quantile(vals, [alpha / 2, 0.5, 1.0 - alpha / 2])
-    return float(q[0]), float(q[1]), float(q[2]), int(len(vals))
-
-
-def s50_mc_ci_from_mvnorm(
-    b,
-    cov,
-    transform="raw",
-    M=MC_DRAWS_TABLE,
-    seed=123,
-    alpha=0.05,
-    enforce_positive_slope=True,
-    slope_eps=1e-10,
-):
-    """Monte Carlo interval for raw-scale midpoint slope s50."""
-    draws, _ = gaussian_parameter_draws(
-        b,
-        cov,
-        M=M,
-        seed=seed,
-        enforce_positive_slope=enforce_positive_slope,
-        slope_eps=slope_eps,
-    )
-
-    if len(draws) == 0:
-        return np.nan, np.nan, np.nan, 0
-
-    vals = np.asarray(
-        [s50_from_b(bb, transform=transform) for bb in draws],
-        float,
-    )
-    vals = vals[np.isfinite(vals)]
-
-    if len(vals) == 0:
-        return np.nan, np.nan, np.nan, 0
-
-    q = np.quantile(vals, [alpha / 2, 0.5, 1.0 - alpha / 2])
-    return float(q[0]), float(q[1]), float(q[2]), int(len(vals))
-
-
-# Historical alias retained.
-def s50_normal_ci_from_mvnorm(
-    b,
-    cov,
-    transform="raw",
-    M=MC_DRAWS_TABLE,
-    seed=123,
-    alpha=0.05,
-    enforce_positive_slope=True,
-    slope_eps=1e-10,
-):
-    return s50_mc_ci_from_mvnorm(
-        b,
-        cov,
-        transform=transform,
-        M=M,
-        seed=seed,
-        alpha=alpha,
-        enforce_positive_slope=enforce_positive_slope,
-        slope_eps=slope_eps,
-    )
-
-
-# ============================================================
-# 8) Bootstrap parameter generators
-# ============================================================
-
-def bootstrap_params_nonparametric(
-    x,
-    y,
-    B=2000,
-    seed=123,
-    l2=0.0,
-    b_start=None,
-):
-    """Ordinary patient-level nonparametric bootstrap."""
-    rng = np.random.default_rng(seed)
-
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-    n = len(y)
-
-    out = []
-    failed = 0
-
-    for _ in range(int(B)):
-        idx = rng.choice(n, size=n, replace=True)
-        xb = x[idx]
-        yb = y[idx]
-
-        if np.unique(yb).size < 2:
-            failed += 1
-            continue
-
-        try:
-            bb = fit_newton(xb, yb, b_start=b_start, l2=l2)
-            if np.all(np.isfinite(bb)):
-                out.append(bb)
-            else:
-                failed += 1
-        except Exception:
-            failed += 1
-
-    arr = np.asarray(out, float) if out else np.empty((0, 2), float)
-    diagnostics = {
-        "attempted": int(B),
-        "successful": int(len(arr)),
-        "failed": int(failed),
-    }
-    return arr, diagnostics
-
-
-
-# Historical alias retained for older notebook cells.
-def bootstrap_params_ordinary(
-    x,
-    y,
-    B=2000,
-    seed=123,
-    l2=0.0,
-    b_start=None,
-):
-    return bootstrap_params_nonparametric(
-        x,
-        y,
-        B=B,
-        seed=seed,
-        l2=l2,
-        b_start=b_start,
-    )
-
-
-def bootstrap_params_stratified(
-    x,
-    y,
-    B=2000,
-    seed=123,
-    l2=0.0,
-    b_start=None,
-):
-    """Class-stratified nonparametric bootstrap preserving class counts."""
-    rng = np.random.default_rng(seed)
-
-    x = np.asarray(x, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-
-    x0 = x[y == 0]
-    x1 = x[y == 1]
-    n0 = len(x0)
-    n1 = len(x1)
-
-    if n0 == 0 or n1 == 0:
-        return np.empty((0, 2), float), {
-            "attempted": int(B),
-            "successful": 0,
-            "failed": int(B),
-        }
-
-    out = []
-    failed = 0
-
-    for _ in range(int(B)):
-        xb0 = rng.choice(x0, size=n0, replace=True)
-        xb1 = rng.choice(x1, size=n1, replace=True)
-
-        xb = np.concatenate([xb0, xb1])
-        yb = np.concatenate([
-            np.zeros(n0, dtype=int),
-            np.ones(n1, dtype=int),
-        ])
-
-        try:
-            bb = fit_newton(xb, yb, b_start=b_start, l2=l2)
-            if np.all(np.isfinite(bb)):
-                out.append(bb)
-            else:
-                failed += 1
-        except Exception:
-            failed += 1
-
-    arr = np.asarray(out, float) if out else np.empty((0, 2), float)
-    diagnostics = {
-        "attempted": int(B),
-        "successful": int(len(arr)),
-        "failed": int(failed),
-        "success_rate": float(len(arr) / B) if B > 0 else np.nan,
-    }
-    return arr, diagnostics
-
-
-def bootstrap_params_parametric(
-    x,
-    b,
-    B=2000,
-    seed=123,
-    l2=0.0,
-    min_ae=2,
-):
-    """
-    Parametric bootstrap with y* ~ Bernoulli[p_hat(x)].
-
-    ``min_ae`` is retained for compatibility.
-    """
-    rng = np.random.default_rng(seed)
-
-    x = np.asarray(x, float).reshape(-1)
-    b = np.asarray(b, float).reshape(2)
-
-    p = model_p(x, b)
-    n = len(x)
-
-    out = []
-    tries = 0
-    max_tries = max(10 * int(B), 1000)
-
-    while len(out) < int(B) and tries < max_tries:
-        tries += 1
-        yb = rng.binomial(1, p, size=n).astype(int)
-
-        n1 = int(np.sum(yb))
-        n0 = n - n1
-
-        if n1 < int(min_ae) or n0 < 1:
-            continue
-
-        try:
-            bb = fit_newton(x, yb, b_start=b, l2=l2)
-            if np.all(np.isfinite(bb)):
-                out.append(bb)
-        except Exception:
-            pass
-
-    arr = np.asarray(out, float) if out else np.empty((0, 2), float)
-    diagnostics = {
-        "attempted": int(tries),
-        "successful": int(len(arr)),
-        "failed_or_rejected": int(tries - len(arr)),
-        "success_rate": float(len(arr) / tries) if tries > 0 else np.nan,
-    }
-    return arr, diagnostics
-
-
-# ============================================================
-# 9) High-level wrapper for one panel
-# ============================================================
-
-def fit_ci_pack_rawgrid(
-    x_raw,
-    y,
-    transform="raw",
-    xmax_raw=None,
-    grid_n=500,
-    name="",
-    l2=0.0,
-    B=2000,
-    seed=123,
-    min_ae=2,
-    z=1.959963984540054,
-):
-    """
-    Fit one validated logistic model and construct five uncertainty summaries.
-    """
-    x_raw = np.asarray(x_raw, float).reshape(-1)
-    y = np.asarray(y, int).reshape(-1)
-
-    if transform not in ("raw", "log"):
-        raise ValueError("transform must be 'raw' or 'log'")
-
-    x_raw = np.clip(x_raw, 1e-12, None)
-    x_model = x_raw if transform == "raw" else np.log(x_raw)
-
-    # Single validated fitting path.
-    validated = fit_pack(
-        x_model,
-        y,
-        name=name,
-        l2=l2,
-        z=z,
-    )
-    b = validated["b"]
-    cov = validated["cov"]
-    gof = validated["gof"]
-
-    xmin_raw = float(np.min(x_raw))
-    xmax0 = float(np.max(x_raw))
-    xmax_use = xmax0 if xmax_raw is None else max(float(xmax_raw), xmax0)
-
-    x_grid_raw = np.linspace(xmin_raw, xmax_use, int(grid_n))
-    x_grid_raw = np.clip(x_grid_raw, 1e-12, None)
-    x_grid_model = x_grid_raw if transform == "raw" else np.log(x_grid_raw)
-
-    mc_x50_bounds = (
-        float(np.min(x_grid_model)),
-        float(np.max(x_grid_model)),
-    )
-
-    mc_draws, diag_mc = gaussian_parameter_draws(
-        b,
-        cov,
-        M=MC_DRAWS_BANDS,
-        seed=seed + 10,
-        enforce_positive_slope=True,
-        x50_bounds=mc_x50_bounds,
-    )
-
-    pars_np, diag_np = bootstrap_params_nonparametric(
-        x_model,
-        y,
-        B=B,
-        seed=seed + 1,
-        l2=l2,
-        b_start=b,
-    )
-
-    pars_str, diag_str = bootstrap_params_stratified(
-        x_model,
-        y,
-        B=B,
-        seed=seed + 2,
-        l2=l2,
-        b_start=b,
-    )
-
-    pars_pm, diag_pm = bootstrap_params_parametric(
-        x_model,
-        b,
-        B=B,
-        seed=seed + 3,
-        l2=l2,
-        min_ae=min_ae,
-    )
-
-    return {
-        "name": name,
-        "transform": transform,
-        "l2": float(l2),
-        "x_raw": x_raw,
-        "x_model": x_model,
-        "y": y,
-        "x_grid_raw": x_grid_raw,
-        "x_grid_model": x_grid_model,
-        "b": b,
-        "cov": cov,
-        "gof": gof,
-        "LCL": validated["LCL"],
-        "UCL": validated["UCL"],
-        "bands": OrderedDict([
-            ("Wald", ci_band_wald(x_grid_model, b, cov, z=z)),
-            ("MC", bootstrap_band_from_params(x_grid_model, mc_draws)),
-            ("Nonparametric", bootstrap_band_from_params(x_grid_model, pars_np)),
-            ("Stratified", bootstrap_band_from_params(x_grid_model, pars_str)),
-            ("Parametric", bootstrap_band_from_params(x_grid_model, pars_pm)),
-        ]),
-        "pars_mc": mc_draws,
-        "pars_nonparametric": pars_np,
-        "pars_stratified": pars_str,
-        "pars_parametric": pars_pm,
-        # Historical aliases
-        "pars_normal": mc_draws,
-        "pars_nonparam": pars_np,
-        "pars_nonparam_ordinary": pars_np,
-        "pars_nonparam_stratified": pars_str,
-        "bootstrap_diagnostics": OrderedDict([
-            ("MC", diag_mc),
-            ("Nonparametric", diag_np),
-            ("Stratified", diag_str),
-            ("Parametric", diag_pm),
-        ]),
-    }
-
-
-# ============================================================
-# 10) Model-band table with LL / UL
-# ============================================================
-
-def model_ci_table_4methods(
-    P,
-    keys=("FULL-RAW", "TRIM-RAW", "FULL-LOG", "TRIM-LOG"),
-):
-    """Long pointwise model-band table."""
-    import pandas as pd
-
-    rows = []
-
-    for key in keys:
-        pk = P[key]
-        xg_raw = np.asarray(pk["x_grid_raw"], float)
-        xg_mod = np.asarray(pk["x_grid_model"], float)
-        trans = pk.get("transform", "")
-
-        for method in CI_METHODS:
-            if method not in pk["bands"]:
-                continue
-
-            lo, md, hi = pk["bands"][method]
-            lo = np.asarray(lo, float)
-            md = np.asarray(md, float)
-            hi = np.asarray(hi, float)
-
-            for i in range(len(xg_raw)):
-                rows.append({
-                    "Panel": key,
-                    "Method": method,
-                    "transform": trans,
-                    "x_grid_raw": float(xg_raw[i]),
-                    "x_grid_model": float(xg_mod[i]),
-                    "fit": float(md[i]),
-                    "LL": float(lo[i]),
-                    "UL": float(hi[i]),
-                    "width": float(hi[i] - lo[i]),
-                })
-
-    return pd.DataFrame(rows)
-
-
-# ============================================================
-# 11) Parameter, x50, and s50 CI summary table
-# ============================================================
-
-def param_ci_table_4methods(
-    P,
-    keys=("FULL-RAW", "TRIM-RAW", "FULL-LOG", "TRIM-LOG"),
-    z=1.959963984540054,
-    alpha=0.05,
-    include_point_est=True,
-    M_normal=MC_DRAWS_TABLE,
-    seed_normal=123,
-):
-    """
-    Notebook-compatible CI summary using:
-      Wald, MC, Nonparametric, Stratified, Parametric.
-    """
-    import pandas as pd
-
-    def _quantile_ci(values):
-        values = np.asarray(values, float)
-        values = values[np.isfinite(values)]
-
-        if len(values) == 0:
-            return np.nan, np.nan, np.nan
-
-        q = np.quantile(values, [alpha / 2, 0.5, 1.0 - alpha / 2])
-        return float(q[0]), float(q[1]), float(q[2])
-
-    def _to_raw_x50_scalar(value, trans):
-        if not np.isfinite(value):
-            return np.nan
-        return float(np.exp(value)) if trans == "log" else float(value)
-
-    def _width(lo, hi):
-        if np.isfinite(lo) and np.isfinite(hi):
-            return float(hi - lo)
-        return np.nan
-
-    rows = []
-
-    for ik, key in enumerate(keys):
-        pk = P[key]
-
-        b = np.asarray(pk["b"], float).reshape(2)
-        cov = np.asarray(pk["cov"], float).reshape(2, 2)
-        trans = pk.get("transform", "")
-        l2 = float(pk.get("l2", 0.0))
-
-        b0_hat = float(b[0])
-        b1_hat = float(b[1])
-        x50_hat = float(x50(b))
-        suv50_hat = _to_raw_x50_scalar(x50_hat, trans)
-        s50_hat = float(s50_from_b(b, transform=trans))
-
-        def add_row(method, b0_ci, b1_ci, x_ci, suv_ci, s_ci, n_used):
-            row = {
-                "Panel": key,
-                "Method": method,
-
-                "b0_hat": b0_hat,
-                "b0_LCL": b0_ci[0],
-                "b0_UCL": b0_ci[2],
-                "b0_width": _width(b0_ci[0], b0_ci[2]),
-
-                "b1_hat": b1_hat,
-                "b1_LCL": b1_ci[0],
-                "b1_UCL": b1_ci[2],
-                "b1_width": _width(b1_ci[0], b1_ci[2]),
-
-                "x50_hat": x50_hat,
-                "x50_med": x_ci[1],
-                "x50_LCL": x_ci[0],
-                "x50_UCL": x_ci[2],
-                "x50_width": _width(x_ci[0], x_ci[2]),
-
-                "SUV50_hat": suv50_hat,
-                "SUV50_med": suv_ci[1],
-                "SUV50_LCL": suv_ci[0],
-                "SUV50_UCL": suv_ci[2],
-                "SUV50_width": _width(suv_ci[0], suv_ci[2]),
-
-                "s50_hat": s50_hat,
-                "s50_med": s_ci[1],
-                "s50_LCL": s_ci[0],
-                "s50_UCL": s_ci[2],
-                "s50_width": _width(s_ci[0], s_ci[2]),
-
-                "transform": trans,
-                "l2": l2,
-                "B_used": n_used,
-            }
-
-            if not include_point_est:
-                for col in (
-                    "b0_hat",
-                    "b1_hat",
-                    "x50_hat",
-                    "x50_med",
-                    "SUV50_hat",
-                    "SUV50_med",
-                    "s50_hat",
-                    "s50_med",
-                    "transform",
-                    "l2",
-                ):
-                    row.pop(col, None)
-
-            rows.append(row)
-
-        # Wald
-        lcl, ucl = wald_ci(b, cov, z=z)
-        x_l, x_u = x50_wald_ci(b, cov, z=z)
-        suv_l = _to_raw_x50_scalar(x_l, trans)
-        suv_u = _to_raw_x50_scalar(x_u, trans)
-        s_l, s_u = s50_wald_ci_numeric(b, cov, transform=trans, z=z)
-
-        add_row(
-            "Wald",
-            (float(lcl[0]), b0_hat, float(ucl[0])),
-            (float(lcl[1]), b1_hat, float(ucl[1])),
-            (float(x_l), x50_hat, float(x_u)),
-            (float(suv_l), suv50_hat, float(suv_u)),
-            (float(s_l), s50_hat, float(s_u)),
-            np.nan,
-        )
-
-        # Draw-based methods
-        mc_draws = pk.get("pars_mc", pk.get("pars_normal"))
-
-        if mc_draws is None or len(mc_draws) < int(M_normal):
-            mc_draws, _ = gaussian_parameter_draws(
-                b,
-                cov,
-                M=M_normal,
-                seed=seed_normal + 1000 * ik,
-                enforce_positive_slope=True,
-            )
-
-        draw_sets = {
-            "MC": mc_draws,
-            "Nonparametric": pk.get(
-                "pars_nonparametric",
-                pk.get("pars_nonparam", np.empty((0, 2))),
-            ),
-            "Stratified": pk.get(
-                "pars_stratified",
-                pk.get("pars_nonparam_stratified", np.empty((0, 2))),
-            ),
-            "Parametric": pk.get(
-                "pars_parametric",
-                np.empty((0, 2)),
-            ),
-        }
-
-        for method in ("MC", "Nonparametric", "Stratified", "Parametric"):
-            pars = np.asarray(draw_sets[method], float)
-
-            if pars.ndim != 2 or len(pars) == 0:
-                nan3 = (np.nan, np.nan, np.nan)
-                add_row(method, nan3, nan3, nan3, nan3, nan3, 0)
-                continue
-
-            b0_ci = _quantile_ci(pars[:, 0])
-            b1_ci = _quantile_ci(pars[:, 1])
-
-            xvals = np.asarray([x50(bb) for bb in pars], float)
-            suvvals = np.asarray(
-                [_to_raw_x50_scalar(v, trans) for v in xvals],
-                float,
-            )
-            svals = np.asarray(
-                [s50_from_b(bb, transform=trans) for bb in pars],
-                float,
-            )
-
-            add_row(
-                method,
-                b0_ci,
-                b1_ci,
-                _quantile_ci(xvals),
-                _quantile_ci(suvvals),
-                _quantile_ci(svals),
-                int(len(pars)),
-            )
-
-    return pd.DataFrame(rows)
-
-
-def combined_x50_model_bounds_table(
-    P,
-    keys=("FULL-RAW", "TRIM-RAW", "FULL-LOG", "TRIM-LOG"),
-    z=1.959963984540054,
-    alpha=0.05,
-    M_normal=MC_DRAWS_TABLE,
-    seed_normal=123,
-):
-    """Combine characteristic intervals with global curve-band summaries."""
-    import pandas as pd
-
-    param_df = param_ci_table_4methods(
-        P,
-        keys=keys,
-        z=z,
-        alpha=alpha,
-        include_point_est=True,
-        M_normal=M_normal,
-        seed_normal=seed_normal,
-    ).copy()
-
-    model_df = model_ci_table_4methods(P, keys=keys).copy()
-
-    global_df = (
-        model_df
-        .groupby(["Panel", "Method"], as_index=False)
-        .agg(
-            global_LL=("LL", "min"),
-            global_UL=("UL", "max"),
-            fit_min=("fit", "min"),
-            fit_max=("fit", "max"),
-            mean_width=("width", "mean"),
-            max_width=("width", "max"),
-        )
-    )
-
-    global_df["global_width"] = global_df["global_UL"] - global_df["global_LL"]
-
-    out = pd.merge(
-        param_df,
-        global_df,
-        on=["Panel", "Method"],
-        how="left",
-    )
-
-    preferred = [
-        "Panel", "Method", "transform",
-        "x50_hat", "x50_LCL", "x50_UCL", "x50_width",
-        "SUV50_hat", "SUV50_LCL", "SUV50_UCL", "SUV50_width",
-        "s50_hat", "s50_LCL", "s50_UCL", "s50_width",
-        "global_LL", "global_UL", "global_width",
-        "fit_min", "fit_max", "mean_width", "max_width", "B_used",
-    ]
-
-    cols = [c for c in preferred if c in out.columns] + [
-        c for c in out.columns if c not in preferred
-    ]
-
-    return out[cols]
-
-
-def bootstrap_diagnostics_table(P):
-    """Return success information for MC and bootstrap methods."""
-    import pandas as pd
-
-    rows = []
-    for panel, pk in P.items():
-        for method, diag in pk.get("bootstrap_diagnostics", {}).items():
-            rows.append({"Panel": panel, "Method": method, **diag})
-
-    return pd.DataFrame(rows)
-
-
-# ============================================================
-# 12) CI figure
-# ============================================================
-def plot_ci_four_panels(P):
-    import numpy as np
-    import matplotlib.pyplot as plt
-    import matplotlib.lines as mlines
-
-    plt.style.use("default")
-
-    COL_NC = "#4c9ed9"
-    COL_AE = "#f28e2b"
-    COL_FIT = "#000000"
-
-    # Methods displayed in the main figure
-    FIGURE_METHODS = [
-        "Wald",
-        "MC",
-        "Nonparametric",
-        "Parametric",
-    ]
-
-    styles = OrderedDict([
-        ("Wald",          ("#2ca02c", "-.", 0.12)),
-        ("MC",            ("#d62728", "--", 0.14)),
-        ("Nonparametric", ("#1f77b4", ":", 0.16)),
-        ("Parametric",    ("#17becf", (0, (6, 2)), 0.14)),
-    ])
-
-    panel_order = [
-        "FULL-LOG",
-        "FULL-RAW",
-        "TRIM-LOG",
-        "TRIM-RAW",
-    ]
-    panel_letters = ["A", "B", "C", "D"]
-
-    fig, axs = plt.subplots(
-        2,
-        2,
-        figsize=(15, 10),
-        dpi=180,
-        sharex="col",
-        sharey=True,
-    )
-
-    for ax, key, letter in zip(
-        axs.flat,
-        panel_order,
-        panel_letters,
-    ):
-        pk = P[key]
-
-        x_raw = np.asarray(pk["x_raw"], float)
-        y = np.asarray(pk["y"], int)
-        xg_raw = np.asarray(pk["x_grid_raw"], float)
-        xg_model = np.asarray(pk["x_grid_model"], float)
-        transform = pk["transform"]
-
-        if transform == "raw":
-            xs = x_raw
-            xg = xg_raw
-        else:
-            xs = np.log(x_raw)
-            xg = xg_model
-
-        rng = np.random.default_rng(123 + ord(letter))
-        jit = (rng.random(len(y)) - 0.5) * 0.04
-
-        ax.scatter(
-            xs[y == 0],
-            (y + jit)[y == 0],
-            s=22,
-            alpha=0.45,
-            color=COL_NC,
-            edgecolors="none",
-            zorder=5,
-        )
-
-        ax.scatter(
-            xs[y == 1],
-            (y + jit)[y == 1],
-            s=24,
-            alpha=0.85,
-            color=COL_AE,
-            edgecolors="none",
-            zorder=5,
-        )
-
-        for method in FIGURE_METHODS:
-            if method not in pk["bands"]:
-                continue
-
-            lo, _, hi = pk["bands"][method]
-            color, linestyle, fill_alpha = styles[method]
-
-            ax.fill_between(
-                xg,
-                lo,
-                hi,
-                color=color,
-                alpha=fill_alpha,
-                zorder=1,
-            )
-
-            ax.plot(
-                xg,
-                lo,
-                color=color,
-                linestyle=linestyle,
-                lw=1.6,
-                zorder=2,
-            )
-
-            ax.plot(
-                xg,
-                hi,
-                color=color,
-                linestyle=linestyle,
-                lw=1.6,
-                zorder=2,
-            )
-
-        fit_curve = model_p(
-            xg_model,
-            pk["b"],
-        )
-
-        ax.plot(
-            xg,
-            fit_curve,
-            color=COL_FIT,
-            lw=2.5,
-            zorder=6,
-        )
-
-        ax.text(
-            0.03,
-            0.95,
-            letter,
-            transform=ax.transAxes,
-            fontsize=15,
-            ha="left",
-            va="top",
-        )
-
-        ax.set_ylim(-0.05, 1.05)
-        ax.grid(False)
-
-        ax.tick_params(
-            axis="both",
-            which="major",
-            labelsize=11,
-            length=4,
-            width=0.8,
-            direction="out",
-        )
-
-    axs[0, 0].set_ylabel(
-        r"$\mathrm{P(AE \mid X = x)}$",
-        fontsize=13,
-    )
-    axs[1, 0].set_ylabel(
-        r"$\mathrm{P(AE \mid X = x)}$",
-        fontsize=13,
-    )
-
-    axs[1, 0].set_xlabel(
-        r"$\log(\mathrm{X})$",
-        fontsize=13,
-    )
-    axs[1, 1].set_xlabel(
-        r"$\mathrm{X}$",
-        fontsize=13,
-    )
-
-    for ax in axs[0, :]:
-        ax.tick_params(
-            axis="x",
-            which="both",
-            labelbottom=False,
-        )
-
-    for ax in axs[:, 1]:
-        ax.tick_params(
-            axis="y",
-            which="both",
-            labelleft=False,
-        )
-
-    labels = {
-        "Wald": "CI: Wald 95%",
-        "MC": "CI: MC propagation 95%",
-        "Nonparametric": "CI: nonparametric bootstrap 95%",
-        "Parametric": "CI: parametric bootstrap 95%",
-    }
-
-    handles = [
-        mlines.Line2D(
-            [],
-            [],
-            marker="o",
-            color=COL_NC,
-            linestyle="None",
-            markersize=8,
-            label="data: NC",
-        ),
-        mlines.Line2D(
-            [],
-            [],
-            marker="o",
-            color=COL_AE,
-            linestyle="None",
-            markersize=8,
-            label="data: AE",
-        ),
-        mlines.Line2D(
-            [],
-            [],
-            color=COL_FIT,
-            lw=2.5,
-            label="fit",
-        ),
-    ]
-
-    for method in FIGURE_METHODS:
-        color, linestyle, _ = styles[method]
-
-        handles.append(
-            mlines.Line2D(
-                [],
-                [],
-                color=color,
-                lw=2,
-                linestyle=linestyle,
-                label=labels[method],
-            )
-        )
-
-    leg = axs[1, 1].legend(
-        handles=handles,
-        loc="lower right",
-        bbox_to_anchor=(0.98, 0.04),
-        fontsize=8.5,
-        frameon=True,
-    )
-
-    leg.get_frame().set_facecolor("white")
-    leg.get_frame().set_edgecolor("#bdbdbd")
-    leg.get_frame().set_linewidth(0.8)
-
-    fig.subplots_adjust(
-        left=0.08,
-        right=0.98,
-        bottom=0.08,
-        top=0.98,
-        wspace=0.06,
-        hspace=0.06,
-    )
-
-    return fig, axs
-# ============================================================
-# ELASTICITY ANALYSIS (x50 and s50)
-# ============================================================
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-
-# ------------------------------------------------------------
-# Core elasticity computation
-# ------------------------------------------------------------
-def elasticity_x50_s50(theta, mode="raw"):
-    """
-    Elasticity for x50 and s50 with respect to theta0 and theta1.
-
-    mode
-    ----
-    'raw' : eta = theta0 + theta1*x
-    'log' : eta = theta0 + theta1*log(x)
-
-    Returns
-    -------
-    dict with:
-        theta0, theta1,
-        x50, s50,
-        E_x50_theta0, E_x50_theta1,
-        E_s50_theta0, E_s50_theta1
-    """
-    theta0, theta1 = map(float, np.asarray(theta).reshape(2))
-
-    if np.abs(theta1) < 1e-12:
-        return dict(
-            theta0=theta0,
-            theta1=theta1,
-            x50=np.nan,
-            s50=np.nan,
-            E_x50_theta0=np.nan,
-            E_x50_theta1=np.nan,
-            E_s50_theta0=np.nan,
-            E_s50_theta1=np.nan,
-        )
-
-    # =========================
-    # RAW MODEL
-    # =========================
-    if mode == "raw":
-        x50 = -theta0 / theta1
-        s50 = theta1 / 4.0
-
-        E_x50_theta0 = 1.0
-        E_x50_theta1 = -1.0
-
-        E_s50_theta0 = 0.0
-        E_s50_theta1 = 1.0
-
-    # =========================
-    # LOG MODEL
-    # =========================
-    elif mode == "log":
-        x50 = float(np.exp(-theta0 / theta1))
-        s50 = theta1 / (4.0 * x50)
-
-        E_x50_theta0 = -theta0 / theta1
-        E_x50_theta1 = theta0 / theta1
-
-        E_s50_theta0 = -E_x50_theta0
-        E_s50_theta1 = 1.0 - E_x50_theta1
-
-    else:
-        raise ValueError("mode must be 'raw' or 'log'")
-
-    return dict(
-        theta0=theta0,
-        theta1=theta1,
-        x50=x50,
-        s50=s50,
-        E_x50_theta0=E_x50_theta0,
-        E_x50_theta1=E_x50_theta1,
-        E_s50_theta0=E_s50_theta0,
-        E_s50_theta1=E_s50_theta1,
-    )
-
-
-# ------------------------------------------------------------
-# Table for 4 panels
-# ------------------------------------------------------------
-def elasticity_table_4panels(P, keys=None, make_plots=True):
-    import pandas as pd
-    import numpy as np
-
-    if keys is None:
-        keys = list(P.keys())
-
-    rows = []
-
-    for key in keys:
-        pk = P[key]
-
-        if "b" not in pk:
-            print(f"[skip] {key}: no fitted parameter key 'b'")
-            continue
-
-        theta = np.asarray(pk["b"], float).reshape(2)
-        transform = pk.get("transform", "raw")
-
-        res = elasticity_x50_s50(theta, mode=transform)
-
-        rows.append({
-            "Panel": key,
-            "transform": transform,
-            **res
-        })
-
-    df = pd.DataFrame(rows)
-
-    if make_plots and len(df) > 0:
-        plot_x50_values(df)
-        plot_s50_values(df)
-        plot_x50_theta1_elasticity(df)
-        plot_s50_theta1_elasticity(df)
-
-    return df
-
-
-# ------------------------------------------------------------
-# Plots
-# ------------------------------------------------------------
-def plot_x50_values(df):
-    fig, ax = plt.subplots(figsize=(7, 4))
-    ax.bar(df["Panel"], df["x50"])
-    ax.set_ylabel("x50")
-    ax.set_title("x50 across panels")
-    plt.xticks(rotation=30)
-    plt.tight_layout()
-    plt.show()
-
-
-def plot_s50_values(df):
-    fig, ax = plt.subplots(figsize=(7, 4))
-    ax.bar(df["Panel"], df["s50"])
-    ax.set_ylabel("s50")
-    ax.set_title("s50 across panels")
-    plt.xticks(rotation=30)
-    plt.tight_layout()
-    plt.show()
-
-
-def plot_x50_theta1_elasticity(df):
-    fig, ax = plt.subplots(figsize=(7, 4))
-    ax.bar(df["Panel"], df["E_x50_theta1"])
-    ax.set_ylabel("Elasticity")
-    ax.set_title("Elasticity of x50 w.r.t. theta1")
-    plt.xticks(rotation=30)
-    plt.tight_layout()
-    plt.show()
-
-
-def plot_s50_theta1_elasticity(df):
-    fig, ax = plt.subplots(figsize=(7, 4))
-    ax.bar(df["Panel"], df["E_s50_theta1"])
-    ax.set_ylabel("Elasticity")
-    ax.set_title("Elasticity of s50 w.r.t. theta1")
-    plt.xticks(rotation=30)
-    plt.tight_layout()
-    plt.show()
-# ------------------------------------------------------------
-# logistic helpers for noise analysis
-# ------------------------------------------------------------
-
-# ============================================================
-# NOISE ANALYSIS FOR LOGISTIC MODEL
-# Correct x50 for RAW and LOG models
-# ============================================================
-
-import os
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib.lines as mlines
-
-
-# ------------------------------------------------------------
-# logistic fit / prediction / x50
-# CONSISTENT WITH MAIN LOGISTIC ANALYSIS
-# ------------------------------------------------------------
-
-def fit_logistic_x(x_raw, y, transform="raw", l2=1e-8):
-    x_raw = np.clip(np.asarray(x_raw, float).ravel(), 1e-12, None)
-    y = np.asarray(y, int).ravel()
-
-    if transform == "raw":
-        x_model = x_raw
-    elif transform == "log":
-        x_model = np.log(x_raw)
-    else:
-        raise ValueError("transform must be 'raw' or 'log'")
-
-    return fit_newton(x_model, y, l2=l2)
-
-
-def predict_curve_x(b, x_grid_raw, transform="raw"):
-    x_grid_raw = np.clip(np.asarray(x_grid_raw, float), 1e-12, None)
-
-    if transform == "raw":
-        x_model = x_grid_raw
-    elif transform == "log":
-        x_model = np.log(x_grid_raw)
-    else:
-        raise ValueError("transform must be 'raw' or 'log'")
-
-    return model_p(x_model, b)
-
-
-def x50_from_b(b, transform="raw"):
-    b = np.asarray(b, float).reshape(2)
-    x50_model = x50(b)
-
-    if not np.isfinite(x50_model):
-        return np.nan
-
-    if transform == "raw":
-        return float(x50_model)
-    elif transform == "log":
-        return float(np.exp(x50_model))
-    else:
-        raise ValueError("transform must be 'raw' or 'log'")
-
-
-def check_noise_x50(pack):
-    b = np.asarray(pack["b_clean"], float).reshape(2)
-    transform = pack["transform"]
-    x50_raw = pack["x50"]
-
-    x50_model = x50_raw if transform == "raw" else np.log(x50_raw)
-    p50 = model_p(np.array([x50_model]), b)[0]
-
-    print(
-        "transform =", transform,
-        "| x50_raw =", x50_raw,
-        "| P(x50) =", p50
-    )
-
-
-# ------------------------------------------------------------
-# noise helpers
-# ------------------------------------------------------------
-
-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)
-
-
-def band_quantiles(curves):
-    C = np.vstack(curves)
-    return np.quantile(C, [0.025, 0.5, 0.975], axis=0)
-
-
-# ------------------------------------------------------------
-# build noise bands
-# ------------------------------------------------------------
-
-def noise_logistic_bands(
-    x_raw,
-    y,
-    transform="raw",
-    sigma_mult=0.129,
-    sigma_add=0.144,
-    x_max=5,
-    grid_n=1000,
-    n_refit=200,
-    n_tta=3000,
-    seed=1234,
-    l2=1e-8,
-):
-    rng = np.random.default_rng(seed)
-
-    x_raw = np.clip(np.asarray(x_raw, float).ravel(), 1e-12, None)
-    y = np.asarray(y).astype(int).ravel()
-
-    xc = np.linspace(1e-12, x_max, grid_n)
-
-    b_clean = fit_logistic_x(x_raw, y, transform=transform, l2=l2)
-    clean = predict_curve_x(b_clean, xc, transform=transform)
-    x50_val = x50_from_b(b_clean, transform=transform)
-
-    curves = []
-    for _ in range(n_refit):
-        xn = add_noise_mult(x_raw, sigma_mult, rng)
-        bn = fit_logistic_x(xn, y, transform=transform, l2=l2)
-        curves.append(predict_curve_x(bn, xc, transform=transform))
-    mult_refit = band_quantiles(curves)
-
-    curves = []
-    for _ in range(n_tta):
-        xn = add_noise_mult(xc, sigma_mult, rng)
-        curves.append(predict_curve_x(b_clean, xn, transform=transform))
-    mult_tta = band_quantiles(curves)
-
-    curves = []
-    for _ in range(n_refit):
-        xn = add_noise_add(x_raw, sigma_add, rng)
-        bn = fit_logistic_x(xn, y, transform=transform, l2=l2)
-        curves.append(predict_curve_x(bn, xc, transform=transform))
-    add_refit = band_quantiles(curves)
-
-    curves = []
-    for _ in range(n_tta):
-        xn = add_noise_add(xc, sigma_add, rng)
-        curves.append(predict_curve_x(b_clean, xn, transform=transform))
-    add_tta = band_quantiles(curves)
-
-    return {
-        "xc": xc,
-        "clean": clean,
-        "x50": x50_val,
-        "b_clean": b_clean,
-        "transform": transform,
-        "l2": float(l2),
-        "mult_refit": mult_refit,
-        "mult_tta": mult_tta,
-        "add_refit": add_refit,
-        "add_tta": add_tta,
-    }
-
-    from scipy.ndimage import gaussian_filter1d
-
-    lo = np.quantile(curves, 0.025, axis=0)
-    md = np.quantile(curves, 0.500, axis=0)
-    hi = np.quantile(curves, 0.975, axis=0)
-
-# smooth boundaries
-    lo = gaussian_filter1d(lo, sigma=8)
-    md = gaussian_filter1d(md, sigma=8)
-    hi = gaussian_filter1d(hi, sigma=8)
-
-    return lo, md, hi
-
-# ------------------------------------------------------------
-# legend
-# ------------------------------------------------------------
-
-def noise_legend_handles():
-    return [
-        mlines.Line2D([], [], marker="o", color="#2b8cbe",
-                      linestyle="None", markersize=7, label="NC data"),
-        mlines.Line2D([], [], marker="o", color="#d7301f",
-                      linestyle="None", markersize=7, label="AE data"),
-        mlines.Line2D([], [], color="black", lw=2.2, label="Initial fit"),
-        mlines.Line2D([], [], color="#1f78b4", lw=6, alpha=0.24,
-                      label="refit band, multiplicative noise"),
-        mlines.Line2D([], [], color="#1f78b4", lw=6, alpha=0.10,
-                      label="fixed-model band, multiplicative noise"),
-        mlines.Line2D([], [], color="#e66101", lw=6, alpha=0.24,
-                      label="refit band, additive noise"),
-        mlines.Line2D([], [], color="#e66101", lw=6, alpha=0.10,
-                      label="fixed-model band, additive noise"),
-        mlines.Line2D([], [], color="#666666", ls="--", lw=1.2,
-                      label=r"$x_{50}$"),
-    ]
-
-
-# ------------------------------------------------------------
-# plot one panel
-# ------------------------------------------------------------
-
-def plot_noise_panel(ax, pack, kind="mult", label="A", X=None, y=None):
-    COL_MULT = "#1f78b4"
-    COL_ADD  = "#e66101"
-    COL_NC   = "#2b8cbe"
-    COL_AE   = "#d7301f"
-
-    xc = pack["xc"]
-    clean = pack["clean"]
-    x50_val = pack["x50"]
-
-    if kind == "mult":
-        refit = pack["mult_refit"]
-        tta = pack["mult_tta"]
-        color = COL_MULT
-    elif kind == "add":
-        refit = pack["add_refit"]
-        tta = pack["add_tta"]
-        color = COL_ADD
-    else:
-        raise ValueError("kind must be 'mult' or 'add'")
-
-    lo_r, _, hi_r = refit
-    lo_t, _, hi_t = tta
-
-    ax.fill_between(xc, lo_t, hi_t, color=color, alpha=0.10, zorder=1)
-    ax.fill_between(xc, lo_r, hi_r, color=color, alpha=0.24, zorder=2)
-
-    ax.plot(xc, lo_r, color=color, lw=1.0, alpha=0.65, zorder=3)
-    ax.plot(xc, hi_r, color=color, lw=1.0, alpha=0.65, zorder=3)
-
-    ax.plot(xc, clean, color="black", lw=2.2, zorder=5)
-    ax.axvline(x50_val, color="#666666", ls="--", lw=1.2, alpha=0.9, zorder=4)
-
-    lo_r_x = np.interp(x50_val, xc, lo_r)
-    hi_r_x = np.interp(x50_val, xc, hi_r)
-    lo_t_x = np.interp(x50_val, xc, lo_t)
-    hi_t_x = np.interp(x50_val, xc, hi_t)
-
-    if X is not None and y is not None:
-        X = np.asarray(X).ravel()
-        y = np.asarray(y).astype(int).ravel()
-
-        ax.scatter(
-            X[y == 0], np.zeros(np.sum(y == 0)),
-            color=COL_NC, s=24, alpha=0.75,
-            edgecolors="none", zorder=7
-        )
-        ax.scatter(
-            X[y == 1], np.ones(np.sum(y == 1)),
-            color=COL_AE, s=24, alpha=0.75,
-            edgecolors="none", zorder=7
-        )
-
-    ax.text(0.03, 0.97, label, transform=ax.transAxes,
-            ha="left", va="top", fontsize=15)
-
-    variant_txt = "FULL" if label in ["A", "B"] else "TRIM"
-    d_ref = hi_r_x - lo_r_x
-    d_tta = hi_t_x - lo_t_x
-
-    info_txt = (
-        f"{variant_txt}\n"
-        f"$x_{{50}}$={x50_val:.2f}\n"
-        f"$\\Delta r$={d_ref:.2f}  $\\Delta t$={d_tta:.2f}"
-    )
-
-    ax.text(
-        0.02, 0.14,
-        info_txt,
-        transform=ax.transAxes,
-        fontsize=10,
-        color="#222",
-        ha="left", va="bottom",
-        bbox=dict(facecolor="white", edgecolor=color,
-                  boxstyle="square,pad=0.25", alpha=0.9)
-    )
-
-    ax.set_xlim(0, xc.max())
-    ax.set_ylim(-0.05, 1.05)
-    ax.grid(alpha=0.25)
-    ax.tick_params(axis="both", labelsize=10)
-
-
-# ------------------------------------------------------------
-# full noise figure
-# ------------------------------------------------------------
-
-def plot_noise_figure(
-    pack_full,
-    pack_trim,
-    X_full,
-    y_full,
-    X_trim,
-    y_trim,
-    figsize=(12, 9),
-    dpi=300,
-):
-    fig, axes = plt.subplots(
-        2, 2,
-        figsize=figsize,
-        dpi=dpi,
-        sharex=True,
-        sharey=True
-    )
-
-    axes = axes.ravel()
-
-    plot_noise_panel(axes[0], pack_full, kind="mult", label="A", X=X_full, y=y_full)
-    plot_noise_panel(axes[1], pack_full, kind="add",  label="B", X=X_full, y=y_full)
-    plot_noise_panel(axes[2], pack_trim, kind="mult", label="C", X=X_trim, y=y_trim)
-    plot_noise_panel(axes[3], pack_trim, kind="add",  label="D", X=X_trim, y=y_trim)
-
-    axes[0].set_ylabel(r"$P(\mathrm{AE}\mid X=x)$", fontsize=12)
-    axes[2].set_ylabel(r"$P(\mathrm{AE}\mid X=x)$", fontsize=12)
-    axes[2].set_xlabel(r"$X$", fontsize=12)
-    axes[3].set_xlabel(r"$X$", fontsize=12)
-
-    handles = noise_legend_handles()
-
-    leg = axes[3].legend(
-        handles=handles,
-        loc="lower right",
-        bbox_to_anchor=(0.97, 0.05),
-        fontsize=9,
-        frameon=True
-    )
-
-    frame = leg.get_frame()
-    frame.set_facecolor("white")
-    frame.set_edgecolor("#bdbdbd")
-    frame.set_linewidth(0.8)
-
-    fig.tight_layout()
-    return fig, axes
-
-
-# ------------------------------------------------------------
-# wrapper
-# ------------------------------------------------------------
-
-def make_noise_figure(
-    X_full,
-    y_full,
-    X_trim,
-    y_trim,
-    transform="raw",
-    sigma_mult=0.129,
-    sigma_add=0.144,
-    x_max=5,
-    grid_n=1000,
-    n_refit=1100,
-    n_tta=10000,
-    seed=1234,
-    l2=1e-8,
-    save_path=None,
-):
-    pack_full = noise_logistic_bands(
-        X_full, y_full,
-        transform=transform,
-        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,
-        l2=l2,
-    )
-
-    pack_trim = noise_logistic_bands(
-        X_trim, y_trim,
-        transform=transform,
-        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,
-        l2=l2,
-    )
-
-    print("FULL n:", len(X_full), "x50:", pack_full["x50"])
-    print("TRIM n:", len(X_trim), "x50:", pack_trim["x50"])
-    check_noise_x50(pack_full)
-    check_noise_x50(pack_trim)
-
-    fig, axes = plot_noise_figure(
-        pack_full, pack_trim,
-        X_full, y_full,
-        X_trim, y_trim,
-        figsize=(12, 9),
-        dpi=300,
-    )
-
-    if save_path is not None:
-        folder = os.path.dirname(save_path)
-        if folder:
-            os.makedirs(folder, exist_ok=True)
-
-        fig.savefig(f"{save_path}.png", dpi=300, bbox_inches="tight")
-        fig.savefig(f"{save_path}.pdf", bbox_inches="tight")
-
-    return fig, axes, pack_full, pack_trim