qra/python/bayesian/Bayesian_Zahra_v1.1.py

# Constrained Bayesian fit (Gamma for NC, Beta-Prime for AE) — no regularization, no r-prior
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy.special import betaln, gammaln
import scipy.io as io
import os

# determine git root
root = os.popen("git rev-parse --show-toplevel").read().strip()

suv   = io.loadmat(os.path.join(root, "data", "suv_percentilesSLOthenUWM.mat"))['lung_SUVperc_COMBINED'][0:58, :, :]
flags = io.loadmat(os.path.join(root, "data", "flags_combined.mat"))['flags'][0:58, 3]   # 0=NC, 1=AE

# Feature X = max SUV_94 per subject; label y = flags
X = np.nanmax(suv[:, :, 94], axis=1).astype(float).ravel()
y = np.asarray(flags, int).ravel()

# Guard for logs
X = np.clip(X, 1e-12, None)
p_emp = float(y.mean())

# Small helpers

def logistic(z):
    z = np.clip(z, -60, 60)
    return 1.0 / (1.0 + np.exp(-z))

def sigmoid(t):
    return 1.0 / (1.0 + np.exp(-t))

def softplus(t):
    t = np.asarray(t, float)
    return np.log1p(np.exp(-np.abs(t))) + np.maximum(t, 0.0)

# dE(x) pieces for log-odds
def dE_ess(x, a, b, s, k, th):
    x = np.asarray(x, float)
    return (a - k) * np.log(x) - (a + b) * np.log1p(x / s) + x / th

def dE_const(a, b, s, k, th):
    return -(a * np.log(s)) - betaln(a, b) + k * np.log(th) + gammaln(k)

def dE_full(x, a, b, s, k, th):
    return dE_ess(x, a, b, s, k, th) + dE_const(a, b, s, k, th)

# Monotonicity cap for theta
def theta_max(a, b, k, s, eps=1e-12):
    A = a - k
    if A <= 0:
        return np.inf
    r = np.sqrt(a + b) - np.sqrt(max(A, eps))
    if r <= 1e-12:
        return np.inf
    return s / (r * r)

# φ = [p_raw, b_raw, s_raw, k_raw, d_raw, u_raw]  (unconstrained)
def unpack_phi_mono(phi):
    p_raw, b_raw, s_raw, k_raw, d_raw, u_raw = phi
    p = sigmoid(p_raw)                       # (0,1)
    b = softplus(b_raw) + 1e-6               # >0
    s = softplus(s_raw) + 1e-6               # >0
    k = softplus(k_raw) + 1e-6               # >0
    delta = softplus(d_raw) + 1e-6           # >0
    a = k + delta                            # enforce a > k
    th_cap = theta_max(a, b, k, s)           # theta cap from monotonicity
    th = th_cap * sigmoid(u_raw)             # 0 < theta <= th_cap
    return p, a, b, s, k, th

#  Prior on p 
TAU = 25.0                                   # shrink toward empirical AE rate
alpha = max(TAU * p_emp, 1e-6)
beta  = max(TAU * (1.0 - p_emp), 1e-6)

prior_r = None
#prior_r = (1.01, 1.01)
#prior_r = (3, 3)

# Objective: negative log-posterior (likelihood + Beta prior on p)
def neg_post_phi_mono(phi, X, y):
    p, a, b, s, k, th = unpack_phi_mono(phi)
    eps = 1e-12

    L  = (np.log(p) - np.log(1 - p)) + dE_full(X, a, b, s, k, th)
    px = logistic(L)
    nll = -np.sum(y * np.log(px + eps) + (1 - y) * np.log(1 - px + eps))

    # Beta(alpha, beta) prior on p → negative log-prior
    npr_p = -((alpha - 1) * np.log(p + eps) + (beta - 1) * np.log(1 - p + eps))

    if prior_r is None: return nll + npr_p
        
    # Prior on r = theta / theta_max (softly avoid boundaries)
    thcap = theta_max(a, b, k, s)
    if np.isfinite(thcap) and thcap > 0:
        r = np.clip(th / thcap, 1e-9, 1 - 1e-9)
        # Negative log Beta prior:  -[(α-1)log r + (β-1)log(1-r)]  (const dropped)
        npr_r = -((prior_r[0] -1) * np.log(r) + (prior_r[1]-1)* np.log(1.0 - r))
    else:
        npr_r = 0.0

    return nll + npr_p + npr_r

# Initialization (stable, simple)
def init_phi(X, y):
    # Gamma(k, theta) MoM for NC group (only need k0 as a safe size proxy)
    X0 = X[y == 0]
    m0 = X0.mean() if X0.size else X.mean()
    v0 = X0.var()  if X0.size else X.var()
    k0 = 2.0 if v0 <= 0 else max((m0**2)/(v0 + 1e-9), 1.5)

    # AE median to seed s0
    X1 = X[y == 1]
    m1 = np.median(X1) if X1.size else np.median(X)

    p0 = np.clip(float(y.mean()), 1e-3, 1 - 1e-3)
    b0, s0 = 1.5, max(m1, 0.5)

    return np.array([
        np.log(p0 / (1 - p0)),             # p_raw
        np.log(np.expm1(b0) + 1e-9),       # b_raw
        np.log(np.expm1(s0) + 1e-9),       # s_raw
        np.log(np.expm1(k0) + 1e-9),       # k_raw
        np.log(np.expm1(1.0) + 1e-9),      # d_raw  (delta)
        -0.2                               # u_raw  (keeps theta a bit below cap initially)
    ], float)

#  Fit wrapper (one retry)
def fit_bayes_mono(X, y, phi_start=None, rng=None):
    if rng is None:
        rng = np.random.default_rng(0)
    if phi_start is None:
        phi_start = init_phi(X, y)

    obj = lambda phi: neg_post_phi_mono(phi, X, y)

    res = optimize.minimize(
        obj, phi_start, method="L-BFGS-B",
        options={"maxiter": 6000, "ftol": 1e-9}
    )
    if not (res.success and np.isfinite(res.fun)):
        phi_try = phi_start + rng.normal(0, 0.2, size=phi_start.shape)
        res = optimize.minimize(
            obj, phi_try, method="L-BFGS-B",
            options={"maxiter": 6000, "ftol": 1e-9}
        )
    return unpack_phi_mono(res.x), res

#  prediction 
def P_with(theta, x):
    p, a, b, s, k, th = theta
    L = (np.log(p) - np.log(1 - p)) + dE_full(x, a, b, s, k, th)
    return logistic(L)

# Run fit + plot
theta_hat, res = fit_bayes_mono(X, y)
print("Optimization success:", res.success, " fval:", float(res.fun))

(p,a,b,s,k,th) = theta_hat
thcap = theta_max(a, b, k, s)
print("theta (p,a,b,s,k,theta):", tuple(float(t) for t in theta_hat), "ratio(th):", th/thcap)

# x-range (cap right end at 10 for readability)
x_lo = max(1e-6, float(X.min()) * 0.8)
x_hi = min(10.0, float(X.max()) * 1.2)
xg   = np.linspace(x_lo, x_hi, 600)
p_curve = P_with(theta_hat, xg)

fig, ax = plt.subplots(figsize=(7.0, 4.6), dpi=140)
ax.plot(xg, p_curve, color="#000000", lw=2.2, label="P(AE|x) (MAP)")

# overlay data with tiny vertical jitter so points don't overlap
rng_plot = np.random.default_rng(999)
jit = (rng_plot.random(len(y)) - 0.5) * 0.06
ax.scatter(X[y==0], (y + jit)[y==0], s=22, alpha=0.55, color="#2ca02c", edgecolors='none', label='NC')
ax.scatter(X[y==1], (y + jit)[y==1], s=26, alpha=0.75, color="#ff7f0e", edgecolors='none', label='AE')

ax.set_ylim(-0.05, 1.05)
ax.set_xlabel('x')
ax.set_ylabel('P(AE | x)')

if prior_r is None:
    ax.set_title('Constrained Bayesian fit (no regularization, prior on p)')
else:
    ax.set_title(f'Constrained Bayesian fit (no regularization, prior on p and prior r{prior_r})')

ax.grid(alpha=0.3)
ax.legend(loc='lower right', frameon=False)
plt.tight_layout()
plt.show()