Adding a new Script by Zahra

src/irae_risk/full_logistic_Zahra

@@ -0,0 +1,2326 @@
+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