Merge branch 'master' of https://git0.fmf.uni-lj.si/horvat/uncertainty_study

Bayesian_Zahra

@@ -0,0 +1,359 @@
+#BetaPrime vs Gamma, hard-mono fit WITH constant + weak regularization
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import optimize
+from scipy.special import betaln, gammaln
+import scipy.io as io
+
+# ---------- Data ----------
+data_path = "../data/"
+
+# === Load data ===
+data_path = "../data/"
+suv = io.loadmat(data_path + "suv_percentilesSLOthenUWM.mat")['lung_SUVperc_COMBINED'][0:58, :, :]
+flags = io.loadmat(data_path + "flags_combined.mat")['flags'][0:58, 3]
+
+X = np.nanmax(suv[:, :, 94], axis=1).reshape(-1)
+X_NC = X[flags == 0]
+X_AE = X[flags == 1]
+
+
+y = np.array([1]*len(X_AE) + [0]*len(X_NC), int) # 1=AE, 0=NC
+X = np.concatenate([X_AE, X_NC], axis=0)
+# sanity checks now that X,y actually exist
+n = len(y); n1 = int(y.sum()); p_emp = n1 / n
+rng = np.random.default_rng(12345)
+
+# ---------- 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)
+
+# --- Eq: logit P(AE|x) = log(p/(1-p)) + dE(x) ---
+# where dE(x) = dE_ess(x) + C(params)
+
+def dE_ess(x, a, b, s, k, th):
+    # Essential (x-dependent) terms:
+    # (a - k) * log(x) - (a + b) * log(1 + x/s) + x/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):
+    # Constant (parameter-only) terms:
+    # C = -a*log(s) - log B(a,b) + k*log(th) + log Γ(k)
+    return -(a*np.log(s)) - betaln(a, b) + k*np.log(th) + gammaln(k)
+
+def dE_full(x, a, b, s, k, th):
+    # Total evidence term: dE(x) = dE_ess(x) + C
+    return dE_ess(x, a, b, s, k, th) + dE_const(a, b, s, k, th)
+
+# ---------- Global 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)
+
+# phi = [p_raw, b_raw, s_raw, k_raw, d_raw, u_raw]
+def unpack_phi_mono(phi):
+    p_raw, b_raw, s_raw, k_raw, d_raw, u_raw = phi
+    
+    # Map raw parameters into valid constrained space:
+    # p = sigmoid(p_raw) ∈ (0,1) → AE prior (class prior)
+    p = sigmoid(p_raw)                     
+    
+    # b = softplus(b_raw) > 0 → Beta–Prime shape parameter
+    b = softplus(b_raw) + 1e-6             
+    
+    # s = softplus(s_raw) > 0 → Beta–Prime scale parameter
+    s = softplus(s_raw) + 1e-6             
+    
+    # k = softplus(k_raw) > 0 → Gamma shape parameter
+    k = softplus(k_raw) + 1e-6             
+    
+    # delta = softplus(d_raw) > 0; a = k + delta > k
+    delta = softplus(d_raw) + 1e-6         
+    a = k + delta                          
+    
+    # θ (theta) is constrained: 0 < θ ≤ θ_max(a,b,k,s)
+    th_cap = theta_max(a, b, k, s)         
+    th = th_cap * sigmoid(u_raw)            # map u_raw ∈ R into (0, th_cap]
+    
+    return p, a, b, s, k, th
+
+# ---------- Priors ----------
+# Beta prior on p centered at empirical rate
+TAU = 25.0   # reduce to ~5 if you want it weaker
+alpha = max(TAU * float(p_emp), 1e-6)
+beta  = max(TAU * (1.0 - float(p_emp)), 1e-6)
+
+# Weak log-normal shrinkage on positive parameters
+def nlog_lognormal(x, mu, sigma, eps=1e-12):
+    # -log LogNormal(x | mu, sigma) up to additive const
+    x = np.maximum(x, eps)
+    lx = np.log(x)
+    return 0.5 * ((lx - mu)/sigma)**2 + lx
+
+# ---------- Objective ----------
+def neg_post_phi_mono_WITH_CONST_REG(phi, X, y):
+    p, a, b, s, k, th = unpack_phi_mono(phi)
+    eps = 1e-12
+
+    # Likelihood with constant included
+    z  = (np.log(p) - np.log(1-p)) + dE_full(X, a, b, s, k, th)
+    px = logistic(z)
+    nll = -np.sum(y*np.log(px + eps) + (1-y)*np.log(1 - px + eps))
+
+    # Prior on p ~ Beta(alpha, beta)
+    npr_p = -((alpha-1)*np.log(p + eps) + (beta-1)*np.log(1 - p + eps))
+
+    # --- Regularization (weak priors) ---
+    # AE median m1: use AE median if present; otherwise overall median.
+    if (y == 1).any():
+        m1 = np.median(X[y == 1])
+    else:
+        m1 = np.median(X)
+
+    reg = 0.0  # total penalty starts at zero
+
+    # 1) Gamma shape k (>0): very weak prior centered at 2 (σ=1.2).
+    reg += nlog_lognormal(k, mu=np.log(2.0), sigma=1.2)
+
+    # 2) Beta-Prime shape b (>0): same weak prior.
+    reg += nlog_lognormal(b, mu=np.log(2.0), sigma=1.2)
+
+    # 3) Beta-Prime scale s (>0): center near AE median (tighter σ=0.5).
+    reg += nlog_lognormal(s, mu=np.log(max(m1, 1e-6)), sigma=0.5)
+
+    # 4) Left-tail gap delta = a - k (>0): center around ~1.5 (σ=0.5)
+    delta = a - k
+    reg += nlog_lognormal(delta, mu=np.log(1.5), sigma=0.5)
+
+    # Keep theta away from the boundary: Beta(3,3) on r = th/th_cap
+    thcap = theta_max(a, b, k, s)
+    if np.isfinite(thcap) and thcap > 0:
+        r = np.clip(th/thcap, 1e-9, 1-1e-9)
+        npr_r = -((3-1)*np.log(r) + (3-1)*np.log(1 - r))
+    else:
+        npr_r = 0.0
+
+    return nll + npr_p + reg + npr_r
+
+# ---------- Initialization ----------
+def init_phi(X, y):
+    # Method-of-moments init for Gamma(k, theta) using NC data (y==0)
+    # ref: https://en.wikipedia.org/wiki/Gamma_distribution#Estimation_of_parameters
+    X0 = X[y==0]
+    m0 = X0.mean() if X0.size else X.mean()  # sample mean
+    v0 = X0.var()  if X0.size else X.var()   # sample variance
+
+    if v0 <= 0:
+        # If variance is degenerate, pick a safe, sane starting point
+        k0, th0 = 2.0, max(m0/2, 0.1)
+    else:
+        # MoM: k = m^2 / v, theta = v / m, with small eps and lower bounds
+        k0  = max((m0**2)/(v0 + 1e-9), 1.5)
+        th0 = max(v0/(m0 + 1e-9), 0.3)
+
+    X1 = X[y==1]
+    m1 = np.median(X1) if X1.size else np.median(X)
+
+
+# Empirical AE rate as a starting prior for p. Clip away from 0/1 so logit is finite.
+# p0 = n_AE / n, but truncated to [1e-3, 1-1e-3] to avoid infinities in log(p/(1-p)).
+    p0 = np.clip(float(y.mean()), 1e-3, 1 - 1e-3)
+    
+    
+    
+# Simple, stable seeds for AE Beta–Prime:
+#   b0 = 1.5   → mild shape; not too spiky, not too flat.
+#   s0 = max(m1, 0.5) → anchor scale near the AE median, but don’t go tiny.
+    b0, s0 = 1.5, max(m1, 0.5)
+
+
+# Pack raw parameters φ for the optimizer.
+# We optimize in an unconstrained space and map with:
+#   p      = sigmoid(p_raw)
+#   b,s,k  = softplus(raw) + 1e-6
+#   a      = k + softplus(delta_raw) + 1e-6
+#   theta  = theta_max * sigmoid(u_raw)
+# To “invert” softplus for the initial guess we use log(expm1(v)) which is the exact inverse
+# of softplus when you define softplus(t) = log(1 + exp(t)). The +1e-9 is just numerical padding.
+    raw = 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),     # delta_raw
+        -0.2                              # u_raw (keeps theta a bit below cap initially)
+    ], float) 
+    return raw
+
+# ---------- Fitting ----------
+def fit_hard_mono_WITH_CONST_REG(X, y, phi_start=None, maxtries=6, jitter=0.3, rng=None):
+    if rng is None:
+        rng = np.random.default_rng(12345)
+    if phi_start is None:
+        phi_start = init_phi(X, y)
+    phi = phi_start.copy()
+    last_err = None
+    for _ in range(maxtries):
+        res = optimize.minimize(
+            neg_post_phi_mono_WITH_CONST_REG, phi, args=(X, y),
+            method="L-BFGS-B",
+            options=dict(maxiter=12000, ftol=1e-10)
+        )
+        if res.success and np.isfinite(res.fun):
+            return unpack_phi_mono(res.x), res
+        last_err = res
+        phi = phi + rng.normal(0, jitter, size=phi.shape)
+    raise RuntimeError(f"Fit failed. Last status: {getattr(last_err, 'message', 'n/a')}")
+
+# ---------- Convenience ----------
+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)
+
+def diag_report(theta, X):
+    p, a, b, s, k, th = theta
+    thcap = theta_max(a, b, k, s)
+    C = dE_const(a, b, s, k, th)
+    logit_p = np.log(p) - np.log(1 - p)
+    A = a - k
+    den = np.sqrt(a + b) - np.sqrt(max(A, 1e-12))
+    xs = np.inf if den <= 1e-12 else s*np.sqrt(max(A,1e-12))/den
+    print({
+        "p": p, "a": a, "b": b, "s": s, "k": k, "theta": th,
+        "theta_max": thcap, "theta/theta_max": (th/thcap if np.isfinite(thcap) else np.nan),
+        "logit(p)": logit_p, "C": C, "x* (bottleneck)": xs
+    })
+
+def plot_s_shape(theta, X, y, rng=None, ax=None, label='P(AE | x)'):
+    if rng is None:
+        rng = np.random.default_rng(0)
+    if ax is None:
+        fig, ax = plt.subplots(figsize=(7, 4.5))
+
+    x_lo = max(1e-6, float(X.min())*0.8)
+    x_hi = float(X.max())*1.2
+    xg = np.linspace(x_lo, x_hi, 600)
+    pg = P_with(theta, xg)
+
+    ax.plot(xg, pg, lw=2, label=label)
+    jit = (rng.random(len(X)) - 0.5) * 0.06
+    y_jit = y + jit
+    ax.scatter(X[y==0], y_jit[y==0], s=22, alpha=0.35, label='NC (y=0)', edgecolors='none')
+    ax.scatter(X[y==1], y_jit[y==1], s=28, alpha=0.60, label='AE (y=1)', edgecolors='none')
+
+    ax.set_ylim(-0.05, 1.05)
+    ax.set_xlim(x_lo, x_hi)
+    ax.set_xlabel('x')
+    ax.set_ylabel('P(AE | x)')
+    ax.set_title('S-shaped P(AE | x) with hard-mono fit (constant included, regularized)')
+    ax.grid(True, alpha=0.3)
+    ax.legend(loc='lower right', frameon=False)
+    return ax
+
+# ---------- Run fit ----------
+theta_hat, res = fit_hard_mono_WITH_CONST_REG(X, y, rng=rng)
+print("Optimization success:", res.success, "fval:", res.fun)
+diag_report(theta_hat, X)
+
+ax = plot_s_shape(theta_hat, X, y, rng=rng)
+plt.show()
+
+
+
+
+#CI Estimation
+
+
+
+# Delta-method 
+import numdifftools as nd  
+
+# wrap scalar objective for numdifftools
+def build_objective(X, y):
+    def f(phi):
+        return neg_post_phi_mono_WITH_CONST_REG(np.asarray(phi, float), X, y)
+    return f
+
+# compute Σ_φ (covariance in phi-space) at MAP using numdifftools.Hessian
+phi_hat = res.x.copy()                    # MAP in raw-phi space 
+f_obj = build_objective(X, y)             # scalar negative log-posterior
+H = nd.Hessian(f_obj, method='central')(phi_hat)
+Sigma_phi = invert_with_eigenfloor(H, floor=1e-6)
+
+
+# 95% Wald band via Delta method
+z = norm.ppf(0.975)  # 1.96 for 95%  ppf stands for percent point function — it’s the inverse CDF
+
+
+def g_px_at_x(x):
+    """Return g(φ) = P(AE | x, φ), so we can get ∇g(φ̂) via numdifftools.Gradient."""
+    #Build a scalar function g(φ) = P(AE | x, φ) for a fixed x.
+    #We return this function so numdifftools.Gradient can compute ∇g(φ̂).
+    def g(phi):
+        p, a, b, s, k, th = unpack_phi_mono(np.asarray(phi, float))
+        L = (np.log(p) - np.log(1 - p)) + dE_full(x, a, b, s, k, th)
+        return logistic(L)
+    return g
+
+# x-grid
+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, 500)
+
+p_hat = np.empty_like(xg)  # point estimate at MAP
+p_lo  = np.empty_like(xg)  # lower 95%
+p_hi  = np.empty_like(xg)  # upper 95%
+
+for i, x in enumerate(xg):
+    gx = g_px_at_x(x)
+    # point estimate at MAP
+    ph = gx(phi_hat)
+    # gradient wrt φ at φ̂ via numdifftools.Gradient
+    grad = nd.Gradient(gx, method='central')(phi_hat)  # shape (d,)
+    # Delta-method variance on probability scale: var ≈ ∇g^T Σ_φ ∇g
+    var = float(grad @ Sigma_phi @ grad)
+    se  = np.sqrt(max(var, 0.0))
+
+    p_hat[i] = ph
+    p_lo[i]  = np.clip(ph - z * se, 0.0, 1.0)
+    p_hi[i]  = np.clip(ph + z * se, 0.0, 1.0)
+
+# plot
+fig, ax = plt.subplots(figsize=(7,4.5 ))
+
+# curve + band (sharp colors)
+ax.plot(xg, p_hat, color="#000000", lw=2.2, label='P(AE|x) at MAP')
+ax.fill_between(xg, p_lo, p_hi, facecolor="#1f77b4", alpha=0.18, label='95% Wald band (Delta)')
+ax.plot(xg, p_lo, color="#1f77b4")
+ax.plot(xg, p_hi, color="#1f77b4")
+
+# overlay data with tiny vertical jitter
+rng_plot = np.random.default_rng(999)
+jit = (rng_plot.random(len(X)) - 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)')
+ax.set_title('Delta–method band using numdifftools Hessian/Gradient')
+ax.grid(alpha=0.3)
+ax.legend(loc='lower right')
+plt.show()
+
+

Task

@@ -0,0 +1,26 @@
+Meeting on 28 June 2025 with Martin
+
+
+Zahra's Task:
+
+1- Check the Scatter plot for parametrs in Bayesian to see if paremters are distributed around MLE to prove asymptotic assumption. 
+2- IF -1/Hessian equalls to cov of MLE.
+3- Explain the code Properly. 
+
+
+
+Meeting on 04/07/2025
+Zahra Should continue Bayesian
+Readig relevant papers
+Plying with constraints in gamma distribution
+Write the notes until now
+
+
+Short Report of this week:
+1- I applied diffrent clssifiers, SVM, KNN, Random Forest, Dicision Tree.
+2- I tried all possible combined Distributions.
+3- 4 out of them seem to be relaible.
+4- I am trying these 4 distributions.
+5- I am reading the relevant papers. 
+
+Meeting on .....

notes/20250606_minutes.txt

@@ -3,8 +3,22 @@ Date: 8.june 2025
 Meeting with Zahra
 
 Dates:
+  28.6 - Meeting about results and Presentation
   30.6. - Zahra has presentation
-  5.7.  - Zahra visits Iran for 2-3 weeks
+  04.7 - Meeting about Bayesian Model
+  11.07- Meeting on Zoom about distributions
+  23.07- Meeting in Person, Zahra and Katja
+  5.08- Meeting on zoom
+  12.8- Meeting in Person
+  15.08- Meeting on Zoom
+  18.8- Meeting on Zoom
+  26.8- Meeting in Person 
+  05/09-Meeting in Person
+  14/09-17/09: Meetin on Zoon
+  21/09-Meeting on Zoom
+  26/09-Meetin in Person
+  
+  
 
 Asking about accuracy measures:
   * for fit we have goodness of fit measures:
@@ -56,6 +70,9 @@ Ref:
    - https://stats.stackexchange.com/questions/354709/sklearn-metrics-accuracy-score-vs-logisticregression-score
    - https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html
    - https://scikit-learn.org/stable/modules/model_evaluation.html#roc-metrics
+   - https://en.wikipedia.org/wiki/Method_of_moments_(statistics)
+   - https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation
+   - https://en.wikipedia.org/wiki/Limited-memory_BFGS
 
 Plan:
 
@@ -67,5 +84,72 @@ Plan:
 
 Ideal:
   * Have some preliminary results for bayesian model until Friday
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  
+  -----------------------------------------------------------------------------------
+  
+  PIPELINE
+  
+  
+# condprob_flowchart.py
+
+def print_flowchart():
+    print("\n" + "="*45)
+    print(" Conditional Probability Modeling Flowchart")
+    print("="*45 + "\n")
+
+    print("       ┌────────────────────────────┐")
+    print("       │ Define the Modeling Goal   │")
+    print("       └────────────┬──────────────┘")
+    print("                    │")
+    print("       ┌────────────▼──────────────┐")
+    print("       │  Choose Meaningful Model  │")
+    print("       │   and Set Assumptions     │")
+    print("       └────────────┬──────────────┘")
+    print("                    │")
+    print("       ┌────────────▼──────────────┐")
+    print("       │     Fit Model to Data     │")
+    print("       └────────────┬──────────────┘")
+    print("                    │")
+    print("       ┌────────────▼──────────────┐")
+    print("       │  Evaluate Model & Loss    │")
+    print("       └────────────┬──────────────┘")
+    print("                    │")
+    print("       ┌────────────▼──────────────┐")
+    print("       │ Quantify Uncertainty (UQ) │")
+    print("       └────────────┬──────────────┘")
+    print("                    │")
+    print("       ┌────────────▼──────────────┐")
+    print("       │  Interpret Results & UQ   │")
+    print("       └────────────┬──────────────┘")
+    print("                    │")
+    print("       ┌────────────▼──────────────┐")
+    print("       │ Apply in Decision Context │")
+    print("       └────────────────────────────┘\n")
+
+    print("  ✔ Each block corresponds to one function or module")
+    print("  ✔ Update each step with your methods and data later")
+    print("  ✔ Plug in visual tools or UQ techniques where needed")
+
+if __name__ == "__main__":
+    print_flowchart()
+
 
 

python/bayesian/Bayesian_Zahra.py

@@ -0,0 +1,176 @@
+# 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
+
+data_path = "../data/"
+suv   = io.loadmat(data_path + "suv_percentilesSLOthenUWM.mat")['lung_SUVperc_COMBINED'][0:58, :, :]
+flags = io.loadmat(data_path + "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 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 = -((0.5) * np.log(r) + (0.5) * np.log(1.0 - r))
+else:
+    npr_r = 0.0
+
+#  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)
+
+# 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))
+
+    return nll + npr_p
+
+# 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))
+print("theta (p,a,b,s,k,theta):", tuple(float(t) for t in theta_hat))
+
+# 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)')
+ax.set_title('Constrained Bayesian fit (no regularization  prior on p and r:alpha and beta=1.5')
+ax.grid(alpha=0.3)
+ax.legend(loc='lower right', frameon=False)
+plt.tight_layout()
+plt.show()

python/bayesian/Bayesian_Zahra_v1.1.py

@@ -0,0 +1,190 @@
+# 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()

python/bayesian/bayesian.py

@@ -0,0 +1,473 @@
+import numpy as np
+import scipy
+
+"""
+    Setting up scipy distribution to be used in Bayesian model
+
+    Input:
+        distr_str: string of scipy distribution
+        n_pars: int, number of parameters
+        parse_pars: function(pars) -> dict
+    
+    Return: 
+        (log_pdf, distr_sample)
+"""
+def setup_scipy_distr(distr_str, n_pars, parse_pars):
+    
+
+    # log of pdf for all choice
+    def log_pdfs(x, pars, choice = 2):
+        logpdf = eval(distr_str).logpdf
+        if choice in [0,1]: return logpdf(x, **parse_pars(pars))
+        return (logpdf(x, **parse_pars(pars)), logpdf(x, **parse_pars(pars[n_pars:])))
+
+    # sample w.r.t. distribution
+    def distr_sample(rng, pars, n, choice = 2):    
+        rvs = eval(distr_str).rvs
+        if choice in [0,1]: return rvs( **parse_pars(pars), size = n, random_state = rng)
+        return np.concatenate((rvs(**parse_pars(pars), size = n[0], random_state = rng), 
+                            rvs(**parse_pars(pars[n_pars:]), size = n[1], random_state = rng)))
+    
+    return log_pdfs, distr_sample
+
+"""
+    Bayesian model describing conditional probability 
+
+        Prob(X|Y = 1) 
+            = Prob(Y=1)p(X|Y=1)/(Prob(Y=0) p(X|Y=0) + Prob(Y=1) p(X|Y=1))
+            = 1 /(1 + O p(X|Y=0)/p(X|Y=1)) 
+            = 1/(1 + exp(-F))
+
+    where F is the decision function
+
+        F =  log(p(X|Y=1)) - log(p(X|Y=0)) + log(O)
+    
+    and O are the odds
+
+        O = Prob(Y=1)/Prob(Y=0)
+"""
+class BayesianModelRegression:
+    
+    """
+        Constructor
+
+        Input:
+            odds: float, ratio Prob(Y=0)/Prob(Y=1)
+            log_pdfs: function (x, pars, choice = 2) 
+                        match choice:
+                            case 0: return log_pdf0
+                            case 1: return log_pdf1
+                            case _: return (log_pdf0, log_pdf1)
+            bounds: tuple of bounds, (bounds0, bound1)
+
+            distr_sample: function (rng, x, pars, n, choice):
+                          generate n sampled os points using pdfs(choice, pars) 
+    """   
+    def __init__(self, odds, log_pdfs, bounds, distr_sample = None):
+        
+        self.odds = odds
+        self.log_odds = np.log(odds)
+        self.log_pdfs = log_pdfs
+        self.bounds = bounds
+        self.distr_sample = distr_sample
+
+    """
+        Calculate decision function:
+
+             decision = log(p(X|Y=1)) - log(p(X|Y=0)) + log(odds)
+
+        Input:
+            x: float or array of floats
+            pars: parameters for log_pdfs
+
+        Return:
+            float or array of float
+    """
+    def decision(self, x, pars):
+
+        # log of pdf for each group
+        lf0, lf1 = self.log_pdfs(x, pars)
+        
+        return lf1 - lf0 + self.log_odds
+    
+    """
+        Calculate model of the conditional probability Prob(X|Y = 1)
+        
+        Input:
+            x: float or array of floats
+            pars: parameters for log_pdfs
+    """    
+    def model(self, x, pars):
+    
+        # decision function
+        F = self.decision(x, pars)
+
+        # calculating model
+        return 1/(1 + np.exp(-F))
+    
+    """
+        Negative Log Likelihood function:
+
+            neg. log likelihood = -sum_i log(Prob(X = x_i, Y = y_i))
+        
+        where
+
+            Prob(X, Y = 1) = 1/(1 + exp(-F))
+            Prob(X, Y = 0) = 1 -  Prob(X, Y = 1) = 1/(1 + exp(+F))
+        
+        with
+
+            log Prob(X, Y = y) = -log(1 + exp(-S(y) F))
+            S(y) =  [ +1 : y = 1
+                    [ -1 : y = 0
+
+        Input:
+            x: array of floats
+            y: array of ints in {0,1}
+            pars: array of floats, model parameters
+        
+        Return:
+            float: negative log likelihood
+
+    """
+    def nllf(self, x, y, pars):
+        
+        # signs for the groups:
+        # group 1 has + sign and group 0 has - sign
+        S = 2.0*y - 1
+
+        # decision function
+        F = self.decision(x, pars)
+        
+        return np.sum(np.log(1 + np.exp(-S*F)))
+    
+
+    """
+        MLE fitting of a distribution,  given by log_pdf, to data x associated 
+        to the group 0 or 1 by maximizing 
+
+            loglikehood_{single group} = sum_i log_pdf(x | pars)
+
+        Input:
+            x: array of floats
+            choice: int in {0,1}, selecting the group
+            method: string in ["local", "diff_evol", "anneal"]
+            seed: int, seed of the random generator
+        
+        Return:
+            {"pars": pars_MLE, "cost": NLLF at pars_MLE}
+    """
+    
+    def fit_distr(self, x, choice, method = "local", seed = 1977):
+        fname = "fit_distr"
+
+        cost = lambda pars: -np.sum(self.log_pdfs(x, pars, choice))
+        bnds = self.bounds[choice]
+
+        match method:
+            case "local":
+                # random parameters from boundaries
+                pars0 = np.random.default_rng(seed).uniform(*zip(*bnds))
+                # use local optimizer
+                res = scipy.optimize.minimize(cost, pars0, bounds = bnds, method = "L-BFGS-B")
+            case "diff_evol":
+                res = scipy.optimize.differential_evolution(cost, bounds = bnds)
+            case "annel":
+                res = scipy.optimize.dual_annealing(cost, bounds = bnds)
+            case _:
+                assert False, f"{fname}::this method does not exist"
+    
+        return {"pars": res.x, "success": res.success, "cost": res.fun} 
+
+
+    """
+        MLE fitting of the Bayesian model:
+
+            pars_MLE = argmin_pars NLLF(pars| x, y)
+
+        Input:
+            x: array of floats
+            y: array of ints in {0,1}
+            method: string in ["local", "diff_evol", "anneal"], optimizer
+        
+        Return:
+            {"pars": pars_MLE, "cost": NLLF at pars_MLE}
+    """
+    def fit(self, x, y, pars0 = None, method = "local"):
+        fname = "fit"
+
+        # defined nllf as function of parameters, data is already included
+        cost = lambda pars: self.nllf(x, y, pars)
+
+        # joint bounds of two groups
+        bnds = np.concatenate(self.bounds)
+        
+        match method:
+            case "local":
+                # estimate initial guess of parameters (for local method)
+                if pars0 is None:
+                    get_pars = lambda choice: self.fit_distr(x[y == choice], choice)["pars"]
+                    pars0 = np.r_[get_pars(0), get_pars(1)]
+          
+                # optimize using local optimizer
+                res = scipy.optimize.minimize(cost, pars0, bounds = bnds, method = "L-BFGS-B")
+            case "diff_evol":
+                res = scipy.optimize.differential_evolution(cost, bounds = bnds)
+            case "anneal":
+                res = scipy.optimize.dual_annealing(cost, bounds = bnds)
+            case _:
+                assert False, f"{fname}::this method does not exist"
+
+        return {"pars": res.x, "success": res.success, "cost": res.fun} 
+
+    """
+        Producing goodness of fit measures:
+    
+            LLF = log_likelihood function
+            AIC = Akaike information criterion
+            BIC = Bayesian information criterion
+            
+        Input:
+            x: array of floats
+            y: array of ints in {0,1}
+            pars: array of floats,  model parameters
+            thresh: float, default 0.5, threshold value for classification
+
+        Return:
+            {"n": n, "k":k, "dof":n-k, 
+            "LLF": log_likelihood, 
+            "AIC": AIC, 
+            "BIC": BIC, 
+            "A": classification accuracy (threshold values = 0.5 prob)}    
+    """
+    def goodness_of_fit(self, x, y, pars, thresh = 0.5):
+
+        # model probabilities
+        p = self.model(x, pars)
+
+        # log likelihood
+        llf = -self.nllf(x, y, pars)
+        
+        # information criteria
+        k, n = len(pars), len(x)
+        AIC = 2*k - 2*llf
+        BIC = k*np.log(n) - 2*llf
+
+        # chi2
+        dof = n - k
+        r = (y - p)/np.sqrt(p*(1-p))
+        chi2 = np.sum(r**2)
+        p_val = scipy.stats.chi2.sf(chi2, dof)
+
+        # using model as classifier
+        matches = y == np.heaviside(p - thresh, 1)
+
+        return {"LLF": llf, "AIC": AIC, "BIC": BIC, 
+                "A" : np.count_nonzero (matches)/n,
+                "chi2": chi2, "p-value(chi2)": p_val,  # not very useful
+                "n": n, "k": k, "dof": dof}
+
+    """
+        Generate m parameters via non-parametric bootstrapping with minimal constraint
+
+            bootstrapped sampled = (xb, yb) sampled with replacement from (x,y)
+
+        Input:
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+            method: string in ["local", "diff_evol", "anneal"], optimizer
+        
+        Return:
+            array of m x len(pars) floats
+    """
+    def get_nonparam_boots_pars(self, x, y, m, seed = 1, method = "local"):
+        fname = "get_nonparam_boots_pars"
+
+        rng = np.random.default_rng(seed)
+
+        # discussing original data 
+        res = self.fit(x, y, method = method)
+        assert res["success"], f"{fname}::fitting original data failed."
+
+        # generate bootstrapped parameters
+        pars = res["pars"] 
+        lst = [pars]
+        n = len(x)
+        
+        while True:
+            
+            # sampling with replacement with restrictions
+            idx = rng.choice(n, n)
+            if np.sum(y[idx]) in [0, n]: continue
+            
+            res = self.fit(x[idx], y[idx], pars0 = pars, method = method)
+            if not res["success"]: continue
+
+            lst.append(res["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)
+
+
+    """
+        Generate m parameters via non-parametric stratified bootstrapping:
+
+            bootstrapped sampled = (xb, yb) sampled with replacement from (x,y) for each groups separately
+
+        meaning         
+          
+            xb = (sampled with replacement from x0, sampled with replacement from x1)
+            yb = (0 ... 0, 1 ... 1)
+                    n0        n1
+        
+        Note samples from each group in yb is constant and same as in y.
+
+        Input:
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+            method: string in ["local", "diff_evol", "anneal"], optimizer
+        
+        Return:
+            array of m x len(pars) floats
+    """
+    def get_nonparam_strat_boots_pars(self, x, y, m, seed = 1, method = "local"):
+        fname = "get_nonparam_strat_boots_pars"
+
+        rng = np.random.default_rng(seed)
+
+        # discussing original data 
+        res = self.fit(x, y, method = method)
+        assert res["success"], f"{fname}::fitting original data failed."
+
+        # separate data of both groups
+        xs = [x[y == i] for i in range(2)]
+        ns = [len(e) for e in xs]
+        
+        # common vector states
+        yb = np.concatenate([np.full(ns[i], i) for i in range(2)])
+        
+        # generate bootstrapped parameters
+        pars = res["pars"]   
+        lst = [pars]
+        while True:
+            # stratified sampling with replacement
+            xb = np.concatenate([rng.choice(xs[i], ns[i]) for i in range(2)])
+        
+            # do fitting
+            res = self.fit(xb, yb, pars0 = pars, method = method)
+            if not res["success"]: continue
+
+            lst.append(res["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)
+    
+    """
+        Generate m parameters via parametric bootstrapping:
+            
+            bootstrapped sample = (x, yb)  yb ~ B(ymodel)
+        
+        Input:
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+            method: 
+        
+        Return:
+            array of m x len(pars) floats
+    """
+    def get_param_boots_pars(self, x, y, m, seed = 1, method = "local"):
+        fname = "get_param_boots_pars"
+
+        rng = np.random.default_rng(seed)
+        
+        # discussing original data 
+        res = self.fit(x, y, method = method)
+        assert res["success"], f"{fname}::fitting original data failed."
+
+        # calculate predicted conditional probabilities 
+        pars = res["pars"]
+        p = self.model(x, pars)
+
+        # generate bootstrapped parameters
+        lst = [pars]
+        n = len(x)
+
+        while True:
+
+            # Generate new binary outcomes from Bernoulli(p_i)
+            y_sim = rng.binomial(n = 1, p = p)
+
+            if np.sum(y_sim) in [0, n]: continue
+            
+            # fit and get new parameter
+            res = self.fit(x, y_sim, pars0 = pars, method=method)
+            if not res["success"]: continue
+
+            lst.append(res["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)
+    
+    """
+        Generate m parameters via parametric stratified bootstrapping by 
+        sampling x from parametrized distributions associated to individual groups:
+
+            bootstrapped sample = (xb, yb)
+                xb = (sampled from distr for x0, sampled from distr for x1)
+                yb = (0 ... 0, 1 ... 1)
+                        n0        n1
+
+        Input:
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+            method: 
+        
+        Return:
+            array of m x len(pars) floats
+    """
+    def get_param_strat_boots_pars(self, x, y, m, seed = 1, method = "local"):
+        fname = "get_param_strat_boots_pars"
+
+        assert self.distr_sample is not None, f"{fname}::distr_sample is not defined"
+
+        rng = np.random.default_rng(seed)
+        
+        # separate data of both groups
+        xs = [x[y == i] for i in range(2)]
+        ns = [len(e) for e in xs]
+        
+        # separate data of both groups
+        pars_g = np.concatenate([self.fit_distr(e, i)["pars"] for i, e in enumerate(xs)])
+        
+        # discussing original data 
+        res = self.fit(x, y, method = method)
+        assert res["success"], f"{fname}::fitting original data failed."
+
+        # common vector states
+        yb = np.concatenate([np.full(ns[i], i) for i in range(2)])
+
+        # generate bootstrapped parameters
+        pars = res["pars"]
+        lst = [pars]
+        
+        while True:
+
+            # Generate new sample of points for each group
+            xb = self.distr_sample(rng, pars_g, ns)
+            
+            # fit and get new parameter
+            res = self.fit(xb, yb, pars0 = pars, method = method)
+            if not res["success"]: continue
+
+            lst.append(res["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)

python/bayesian/data_utils.py

@@ -0,0 +1,59 @@
+import numpy as np
+
+
+"""
+  Extract data fom dictionaries for specific organ
+
+  Input: 
+    organ: string in ["lung", "bowel", "thyroid"]
+    perc: int, percentiles [1, ... , 100]
+    suv_dict: dict containing percentiles of suv
+    flags_dict: dict containing states
+"""
+
+def get_data(organ, perc, suv_dict, flags_dict, nr_patient = 58):
+    fname = "get_data"
+
+    # get concrete data set
+    suv = suv_dict[organ + '_SUVperc_COMBINED'][:nr_patient,:,:]  # suv percentiles
+    
+    # index of percentile
+    perc_idx = perc -1
+
+    # computing max SUV percentile per patient, ignoring nans
+    x = np.nanmax(suv[:,:,perc_idx], axis = 1)
+
+    # determining index in the flags based on organ
+    match organ:
+        case "lung":
+            flags_idx = 3
+        case "bowel":
+            flags_idx = 1
+        case "thyroid":
+            flags_idx = 5
+        case _:
+            assert False, f"{fname}::this organ {organ = } is not supported"
+
+    # getting state of patients: 0 == NC, 1  == AE
+    y = flags_dict['flags'][:nr_patient, flags_idx]  
+
+    return x, y
+
+
+"""
+    Check if a vector lies within the specified bounds for each dimension.
+
+    Parameters:
+    - vector (np.ndarray): 1D array representing the point to check. Shape: (n,)
+    - bounds (np.ndarray): 2D array of shape (n, 2), where each row is (min, max) for a dimension.
+
+    Returns:
+    - bool: True if the vector is within bounds in all dimensions, False otherwise.
+"""
+def within_bounds(vector: np.ndarray, bounds: np.ndarray) -> bool:
+
+    if vector.shape[-1] != bounds.shape[0]:
+        raise ValueError("Dimension mismatch: vector length and bounds rows must be equal.")
+    
+    return np.apply_along_axis(lambda x: np.all((x >= bounds[:, 0]) & (x <= bounds[:, 1])), -1, vector)
+

python/logistic/logit_utils.py

@@ -2,6 +2,8 @@ import numpy as np
 import scipy
 import scipy.stats
 
+from sklearn import linear_model
+
 """
     Model function 
 
@@ -49,20 +51,22 @@ def logit_poly_model(x, b):
     Input:
         lm: instance linear_model.LogisticRegression
         x: array of n floats
-        y: array of n floats
+        y: array of n int in {0,1}
         d: degree of decision function
     
     Return:
         params = [b_0, ..., b_degree], array of r = degree + 1 floats
 """
 
-def logit_poly_fit(lm, x, y, degree = 1):
+def logit_poly_fit(x, y, degree = 1, lm  = None):
    
-   X_feature = np.column_stack([x**i for i in range(1, degree+1)])
+    if lm is None: lm = linear_model.LogisticRegression(penalty = None)
 
-   lm.fit(X_feature, y)
+    X_feature = np.column_stack([x**i for i in range(1, degree+1)])
 
-   return np.r_[lm.intercept_[0], lm.coef_[0,:]]
+    lm.fit(X_feature, y)
+
+    return np.r_[lm.intercept_[0], lm.coef_[0,:]]
 
 """
     Producing goodness of fit measures:
@@ -263,26 +267,28 @@ def logit_poly_model_quantiles_delta(x, probs, mean_b, cov):
     boostrapped sample = (xb, yb)  by sampling with replacement pairs (x_i, y_i) 
                                    with condition that yb can not be just 0 or just 1 
 
-     Input:
-        lm: linear_model.LogisticRegression
+    Input:
         x : array of n floats
         y : array of n int in {0,1}
         m: integer, number of samples
         degree: int,  degree of decision function
         seed : int, seed for the random generator
-    
+        lm: linear_model.LogisticRegression
+  
     Return:
         array of mx(degree + 1)
     
     Return:
         array of mx(degree + 1)
 """
-def get_nonparam_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
+def get_nonparam_boots_pars(x, y, m, degree = 1, seed = 1977, lm = None):
+
+    if lm is None: lm = linear_model.LogisticRegression(penalty = None)
 
     rng = np.random.default_rng(seed)
     
     # fitting original data
-    pars = logit_poly_fit(lm, x, y, degree = degree)
+    pars = logit_poly_fit(x, y, degree = degree, lm = lm)
     
     n = len(x)
 
@@ -294,7 +300,7 @@ def get_nonparam_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
         idx = rng.choice(n, n)
         if np.sum(y[idx]) in [0, n]: continue
 
-        lst.append(logit_poly_fit(lm, x[idx], y[idx], degree=degree))
+        lst.append(logit_poly_fit(x[idx], y[idx], degree=degree, lm = lm))
         if len(lst) == m: break
 
     return np.array(lst)
@@ -309,22 +315,24 @@ def get_nonparam_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
             yb = (0...0, 1...1)
 
     Input: 
-        lm: linear_model.LogisticRegression
         x : array of n floats
         y : array of n int in {0,1}
         m: integer, number of samples
         degree: int,  degree of decision function
         seed : int, seed for the random generator
+        lm: linear_model.LogisticRegression
     
     Return:
         array of mx(degree + 1)
 """
-def get_nonparam_stratified_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
+def get_nonparam_stratified_boots_pars(x, y, m, degree = 1, seed = 1977, lm = None):
+
+    if lm is None: lm = linear_model.LogisticRegression(penalty = None)
 
     rng = np.random.default_rng(seed)
 
     # pars of original data
-    pars = logit_poly_fit(lm, x, y, degree = degree)
+    pars = logit_poly_fit(x, y, degree = degree, lm = lm)
     
     # statistics about groups
     xs = [x[y == i] for i in range(2)]
@@ -339,7 +347,7 @@ def get_nonparam_stratified_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
         # stratified sampling with replacement
         xb = np.concatenate([rng.choice(xs[i], ns[i]) for i in range(2)])
         # do fitting
-        lst.append(logit_poly_fit(lm, xb, yb, degree = degree))
+        lst.append(logit_poly_fit(xb, yb, degree = degree, lm = lm))
 
     return np.array(lst)
 
@@ -351,13 +359,13 @@ def get_nonparam_stratified_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
     where B is Bernoulli distribution
 
     Input: 
-        lm: linear_model.LogisticRegression
         x : array of n floats
         y : array of n int in {0,1}
         m: integer, number of samples
         degree: int,  degree of decision function
         seed : int, seed for the random generator
-
+        lm: linear_model.LogisticRegression
+    
     Return:
         array of mx(degree + 1)
 
@@ -366,10 +374,12 @@ def get_nonparam_stratified_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
       * https://www.scirp.org/journal/paperinformation?paperid=70962
       * https://en.wikipedia.org/wiki/Bernoulli_distribution
 """
-def get_parametric_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
+def get_parametric_boots_pars(x, y, m, degree = 1, seed = 1977, lm = None):
+
+    if lm is None: lm = linear_model.LogisticRegression(penalty = None)
 
     # first discuss original dataset
-    pars = logit_poly_fit(lm, x, y, degree = degree)
+    pars = logit_poly_fit(x, y, degree = degree, lm = lm)
     p = logit_poly_model(x, pars)
 
     rng = np.random.default_rng(seed)
@@ -384,7 +394,7 @@ def get_parametric_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
 
         if np.sum(y_sim) in [0, n]: continue
         
-        lst.append(logit_poly_fit(lm, x, y_sim, degree = degree))
+        lst.append(logit_poly_fit(x, y_sim, degree = degree, lm = lm))
 
         if len(lst) == m: break
 

python/logistic/logit_utils_gen.py

@@ -0,0 +1,769 @@
+import numpy as np
+import scipy
+import scipy.optimize
+
+import mono_cubic2 as mc
+
+def resize_with_const(v, n, val=0):
+    """
+    Resize a 1D vector to a specified length `n`.
+
+    If the input vector `v` is longer than `n`, it is truncated.
+    If it is shorter, it is padded with the constant value `val`.
+    If it is already of length `n`, it is returned unchanged.
+
+    Parameters:
+        v (array-like): Input 1D vector (list or NumPy array).
+        n (int): Target length of the output vector.
+        val (scalar, optional): Value used to pad if `v` is shorter than `n`. Default is 0.
+
+    Returns:
+        np.ndarray: Resized 1D NumPy array of length `n`.
+    """
+    v = np.asarray(v)
+    
+    if len(v) == n: return v
+    
+    if len(v) > n: return v[:n]
+    
+    return np.concatenate([v, np.full(n - len(v), val)])
+
+def safe_exp(x, max_exp = 700):
+    """
+    A numerically robust version of np.exp that avoids overflow by clipping the input.
+    
+    Parameters:
+        x : array_like
+            Input value or array.
+        max_exp : float
+            Maximum allowed exponent value. np.exp(709) ≈ 8.2e307 (close to float64 max).
+    
+    Returns:
+        array_like
+            The exponential of the input with overflow protection.
+    """
+
+    return np.exp(np.clip(x, -max_exp, max_exp))
+
+def safe_expit(x, max_exp = 700): return 1/(1 + safe_exp(-x, max_exp))
+
+"""
+    Fitting data
+
+        {(x_i, y_i) : x_i in R, y_i in {0,1}, i = 0, ..., n-1 }
+    
+    to model function 
+
+        f(x|pars) = 1/(1 + exp(-F(x|beta)))         beta = beta(pars)
+
+    with log odds of polynomial form:
+
+        log(f(x|pars)/(1 - f(x|pars))) = F(x|beta)  beta = beta(pars)
+
+    where F is decision function (aka logit)
+        
+        F(x|beta) = sum_{i=0}^degree beta_i x^i
+                
+    and coefficients
+        
+        beta = [beta_i(pars)]_{i=0}^degrees
+    
+    Conditional probability
+
+        Prob(Y = y_i|x, pars) = 1/(1 + exp(-s(y) F(x|beta(pars))) 
+
+    with
+
+            s(y) = 2*y -1
+"""
+class LogisticPolyRegression:
+
+    """
+        Constructor
+
+        Input:
+            degree: int, degree of polynomial
+            mono: boolean, default False
+            lambda: None or tuple float, L1 and L2 regularization
+
+    """
+    def __init__(self, degree = 1, mono = False, lam = None):
+
+        self.degree = degree
+        self.mono = mono
+        self.big = 200
+        self.small = 1e-8
+        self.lam = lam
+
+        if mono:
+            assert self.degree in [1, 3], f"Degree {self.degree} not supported in mono!"
+
+        self.mono3 = self.mono and (self.degree == 3)
+
+    """
+        Mapping regression parameters pars to coefficients beta 
+            
+            beta = beta(pars)
+        
+        used in decision function:
+
+            F(x|beta) = sum_{i=0}^degree beta_i x^i
+
+        Input:
+            pars
+        
+        Return:
+            beta
+    """
+    def get_beta(self, pars):
+        
+        beta = mc.forward_map(pars) if self.mono3 else pars
+        return np.array(beta)
+    
+    """
+        Mapping beta to regression parameters used in decision function.
+
+        Input:
+            beta: coefficient beta
+        
+        Return:
+            pars
+    """
+    def get_pars(self, beta):
+
+        pars = mc.backward_map(beta) if self.mono3 else beta
+        return np.array(pars)
+
+    """
+        Calculate jacobian between beta and  regression parameters
+            
+            J = d(beta)/d(pars)
+              = [d(beta_i)/d(pars_a)]_{i,a}
+        and
+            
+            H = [d^2 beta_i/(d(pars_a) d(pars_b))]_{i,a,b}
+
+        Input:
+            pars
+            hess: boolean, False
+        
+        Return:
+            J       if hess = True
+            (J, H)  if hess = False
+    """
+    def get_jac_beta(self, pars, hess = False):
+        
+        n = len(pars)
+            
+        J = mc.forward_map_jacobian(pars) if self.mono3 else np.eye(n)
+            
+        if not hess: return J
+        
+        H = mc.forward_map_hessian(pars) if self.mono3 else np.zeros(shape = (n, n, n))
+
+        return (J, H)
+
+
+    """
+        Model function 
+
+            f(x|pars) = 1/(1 + exp(-F(x|beta)))  beta = beta(pars)
+
+        Input:
+            x: scalar value or a array of values
+            pars: array of r = degree + 1 floats, model parameters array of floats
+
+        Return:
+            model function values 
+    """
+    def model(self, x, pars):
+
+        beta = self.get_beta(pars)  # beta
+        X = np.column_stack([x**i for i in range(len(beta))])
+        F = X @ beta              # decision function, X beta
+
+        return safe_expit(F)
+
+    """
+        Calculate negative log-likelihood function
+
+            nllf = -sum_i log(Prob(Y = y_i| x_i, pars))   : neg. log likelihood
+            
+            grad = [d(nllf)/d(pars_a)]_a                  : jacobian
+        
+        where
+
+            Prob(Y = y_i| x_i, pars) = 1/(1 + exp(-s_i F(x_i| beta)))  beta=beta(pars)
+            s_i = 2*y_i  - 1
+     
+        Input:
+            x: array of n floats
+            y: array of n int in {0,1}  
+            pars: array of r = degree + 1 floats, model parameters array of floats
+            jac: boolean, default False, if jacobian is needed
+        
+        Return:
+            nllf           : if jac is false
+            (nllf, grad)   : if jac is true
+    """
+    def nllf(self, x, y, pars, jac = False):
+
+        beta = self.get_beta(pars)
+        X = np.column_stack([x**i for i in range(len(beta))])
+    
+        # signs
+        s = 2.0*y - 1
+        
+        # decision function for conditional probability Prob(Y = y| x)
+        F = s*(X @ beta)
+
+        nllf = np.sum(np.log(1 + safe_exp(-F))) 
+
+        if not jac: return nllf
+        
+        J = self.get_jac_beta(pars)
+
+        grad = -(s*safe_expit(-F)) @ (X @ J)
+
+        return (nllf, grad)
+
+    """
+        Penalty function
+
+         Input:
+            pars: array of r = degree + 1 floats, model parameters array of floats
+            jac: boolean, default False, if jacobian is needed
+        
+        Return:
+            val           : if jac is false
+            (val, grad)   : if jac is true
+    """
+    def penalty(self, pars, jac = False):
+        
+        val = self.lam[0]*np.sum(np.abs(pars)) + self.lam[1]*np.sum(pars**2)
+
+        if jac:
+            grad = self.lam[0]*np.sign(pars) + 2*self.lam[1]*pars
+            return (val, grad)
+           
+        return val
+
+    """
+        Cost function
+    """
+    def cost(self, x, y, pars, jac = False):
+
+        val = self.nllf(x, y, pars, jac)
+        
+        if self.lam is not None:
+            pen = self.penalty(pars, jac)
+            return (val[0] + pen[0], val[1] + pen[1]) if jac else val + pen
+        
+        return val
+           
+    
+    """
+        Estimate parameters.
+
+        Input:
+            x: array of n floats
+            y: array of n int in {0,1}
+        
+        Return:
+            pars0
+    """
+    def get_est_pars(self, x, y):
+
+        L = np.log(2*len(x) + 1)
+
+        if self.mono3:    
+            z = L*np.polyfit(x, 2.0*y - 1, 1)[::-1]
+            beta = resize_with_const(z, self.degree + 1, 1e-8)    
+        else:
+            beta = L*np.polyfit(x, 2.0*y - 1, self.degree)[::-1]
+
+        return self.get_pars(beta)
+    
+    """
+        Performing logistic regression with log odds of polynomial form:
+
+            log(f(x|pars)/(1 - f(x|pars))) = F(x|beta)  beta = beta(pars)
+        
+        and this gives           
+            
+            f(x|pars) = 1/(1 + exp(-F(x|beta)))
+
+        where coefficient beta = [beta_i(pars)]_{i=0}^degree, with decision 
+        function (aka logit)
+
+            F(x|beta) = sum_{i=0}^degree beta_i x^i
+                
+        Input:
+            x: array of n floats
+            y: array of n int in {0, 1}
+        
+        Return:
+            pars
+    """
+    def fit(self, x, y, pars0 = None, method = "local"):
+        
+        bnds = [(-self.big, self.big)]*(self.degree + 1)
+  
+        if method == "local":
+
+            pars0 = self.get_est_pars(x, y)
+            cf = lambda pars: self.cost(x, y, pars, jac = True)
+            res = scipy.optimize.minimize(cf, x0 = pars0, method = 'L-BFGS-B', 
+                                          jac = True, bounds = bnds, tol=1e-12)
+
+        elif method == "diff_evol":
+            
+            cf = lambda pars: self.cost(x, y, pars, jac = False)
+            res = scipy.optimize.differential_evolution(cf, bounds = bnds, 
+                                                        tol = 1e-8, polish = False)
+
+            cf = lambda pars: self.cost(x, y, pars, jac = True)
+            res = scipy.optimize.minimize(cf, x0 = res.x, method = 'L-BFGS-B', 
+                                          jac = True, bounds = bnds, tol=1e-12)
+
+        elif method == "anneal":
+
+            cf = lambda pars: self.cost(x, y, pars, jac = False)
+            res = scipy.optimize.dual_annealing(cf, bounds = bnds)
+        
+        else:
+            assert False, "This method is not supported."
+        
+        return {"pars": res.x, "cost": res.fun, "success": res.success}
+
+    """
+        Producing goodness of fit measures:
+        
+            LLF = log_likelihood function
+            AIC = Akaike information criterion
+            BIC = Bayesian information criterion
+
+        Input:
+            x: array of n floats
+            y: array of n int in {0,1}
+            pars: array of r = degree+1 floats, model parameters
+            thresh: float, default 0.5, threshold value for classification
+
+        Return:
+            {"n": n, "k":k, "dof":n-k, 
+            "LLF": log_likelihood, 
+            "AIC": AIC, 
+            "BIC": BIC, 
+            "A": classification accuracy (threshold values = 0.5 prob)}
+        
+        Ref:
+            https://en.wikipedia.org/wiki/Logistic_regression
+            https://en.wikipedia.org/wiki/Akaike_information_criterion
+            https://www.medicine.mcgill.ca/epidemiology/joseph/courses/epib-621/logfit.pdf
+    """
+    def goodness_of_fit(self, x, y, pars, thresh = 0.5):
+        
+        # model probabilities
+        p = self.model(x, pars)
+
+        # log likelihood
+        llf = -self.nllf(x, y, pars)
+        
+        # information criteria
+        k, n = len(pars), len(x)
+        AIC = 2*k - 2*llf
+        BIC = k*np.log(n) - 2*llf
+
+        # chi2
+        dof = n - k
+        r = (y - p)/np.sqrt(p*(1-p) + self.small)
+        chi2 = np.sum(r**2)
+        p_val = scipy.stats.chi2.sf(chi2, dof)
+
+        # using model as classifier
+        matches = y == np.heaviside(p - thresh, 1)
+
+        # accuracy A
+        A = np.count_nonzero(matches)/n
+
+        return {"LLF": llf, 
+                "AIC": AIC, 
+                "BIC": BIC, 
+                "A" : A, 
+                "chi2": chi2, 
+                "p-value(chi2)": p_val,  # not very useful
+                "n": n, "k": k, "dof": dof}
+
+    """
+        Calculation of asymptotic variance-covariance matrix of parameters pars
+
+           cov_{asymp}[pars] = H^{-1}
+
+        where H is hessian of nllf 
+
+           H = [d^2(nllf)/(d(pars)_a d(pars)_b ]_{a,b}
+        
+        for the logistic regression of the polynomial model:
+        
+            log(f(x)/(1 - f(x))) ~ sum_{i=0}^degree b_i(pars) x^i
+
+        Input:
+            x: array of n floats
+            pars: array of r = degree+1 floats, model parameters
+            
+        Return:
+            array of rxr floats; r = degree + 1
+        
+        Ref:
+            https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients
+            https://goodboychan.github.io/machine_learning/2020/09/14/02-Regularized-likelihood-methods.html
+    """
+    def cov(self, x, y, pars):
+        
+        # coefficients
+        beta = self.get_beta(pars)
+        
+        # design matrix -- add column of 1's at the beginning of your X_train matrix
+        X = np.column_stack([x**i for i in range(len(beta))])
+
+        # Jacobian J = [dbeta_i/dpars_j]_{ij}
+        J, H = self.get_jac_beta(pars, hess = True)
+
+        # signs
+        s = 2.0*y - 1
+        
+        # decision function for conditional probability Prob(Y = y| x)
+        F = s*(X @ beta)
+
+        # probabilities p_i = P(Y=y_i | x_i)
+        p = safe_expit(F)
+        q = 1 - p
+
+        # calculate hessian
+        L = X @ J
+        H = (L.T*(q*p))@L - np.tensordot((s*q)@X, H, axes = ([0], [0]))
+
+        if self.lam is not None: 
+            Hp = H + 2*self.lam[1]*np.eye(len(pars))   # H' = H + lambda id 
+            iHp = np.linalg.inv(Hp)                    # inv(H')
+
+            return iHp@H@iHp
+        
+        # covariance matrix C_params = H^-1
+        return np.linalg.inv(H)
+
+
+    """
+        Calculating quantiles of the model parameters at given probabilities p 
+        for normal distribution of parameters:
+
+            pars ~ N(mean_pars, cov_pars)
+        
+        Input:
+            probs: array of m floats, probabilities
+            mean_pars: array of r = degree+1 floats, mean model parameters
+            cov_pars: array of rxr floats, variance-covariance matrix of parameters
+        
+        Return:
+            array of mxr floats
+    """
+    def get_pars_quantiles_normal(self, probs, mean_pars, cov_pars):
+        
+        # mean and standard variance parameters
+        locs = mean_pars
+        scales = np.sqrt(np.diag(cov_pars))
+
+        # computing quantiles of parameters
+        return locs + np.outer(scipy.stats.norm.ppf(probs), scales)
+
+    """
+        Calculating quantiles of the model values 
+        
+            f(x|pars) = 1/(1 + exp(-F(x|beta)))     beta = beta(pars)
+        
+        with
+
+            F(x|beta) = sum_{i=0}^degree x^i beta_i
+            
+        at given probabilities p and values x assuming 
+        normal distribution of parameters:
+
+            pars ~ N(mean_pars, cov_pars)
+        
+        This distribution is asymptotic MLE distribution of parameters.
+
+        Input:
+            x: array of n float
+            probs: array of m floats, probabilities
+            mean_pars: array of r = degree+1 floats, mean model parameters
+            cov_pars: array of rxr floats, variance-covariance matrix of parameters
+        
+        Return:
+            array of mxn floats
+    """
+
+    def get_model_quantiles_normal(self, x, probs, mean_pars, cov_pars, 
+                                   exact = True, seed = 1977, m = 10**5):
+        
+        mean_beta = self.get_beta(mean_pars)
+        X = np.column_stack([x**i for i in range(len(mean_beta))])
+        
+        if exact and self.mono3:
+
+            # init random generator
+            rng = np.random.default_rng(seed)
+            pars = rng.multivariate_normal(mean_pars, cov_pars, size = m)
+
+            # get betas
+            beta = np.apply_along_axis(self.get_beta, 1, pars)
+            
+            # quantiles of decision function
+            Q = np.quantile(X@beta.T, probs, axis = 1)
+            
+            return np.apply_along_axis(safe_expit, 1, Q)
+
+        # J = d(beta)/d(pars)
+        J = self.get_jac_beta(mean_pars)
+        
+        # transform data
+        S = X@J
+
+        # mean and standard variance of logit (aka log of odds)
+        locs = X@mean_beta
+        scales = np.sqrt(np.diag(S@cov_pars@S.T))
+
+        # computing quantiles of logit
+        Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)
+        
+        # convert logit to expit
+        return safe_expit(Q)
+
+    """
+        Calculating quantiles using delta method of the model values 
+        
+            f(x|pars) = 1/(1 + exp(-F(x|beta))) beta = beta(pars)
+
+        with
+
+            F(x|beta) = sum_{i=0}^degree x^i beta_i
+
+        at given probabilities p and values x assuming 
+        normal distribution of parameters :
+
+            pars ~ N(mean_pars, cov_pars)
+
+        This distribution is asymptotic MLE distribution of parameters.
+        We approximate exact model with linear expansion
+
+            f(x|pars) = p(x|pars_mean) + dp/db(x|pars_mean) (pars - pars_mean)
+        
+        and the last term is normally distributed.
+
+        Input:
+            x: array of n float
+            probs: array of m floats, probabilities
+            mean_pars: array of r = degree+1 floats, mean model parameters
+            cov_pars: array of rxr floats, variance-covariance matrix of parameters 
+        
+        Return:
+            array of mxn floats
+    """
+    def get_model_quantiles_delta(self, x, probs, mean_pars, cov_pars):
+        
+        mean_beta = self.get_beta(mean_pars)
+        X = np.column_stack([x**i for i in range(len(mean_beta))])
+        F = X@mean_beta
+
+        # J = d(beta)/d(pars)
+        J = self.get_jac_beta(mean_pars)
+
+        # S = d(F)/d(pars)
+        S = X@J
+
+        # attributes of normal distribution of model values
+        locs = safe_expit(F)
+        scales = np.sqrt(np.diag(S@cov_pars@S.T))/(4*np.cosh(F/2)**2)
+
+        # computing quantiles of logit
+        Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)
+
+        return np.clip(Q, a_min = 0, a_max = 1)
+
+    """
+        Generate parameters assuming normal distribution.
+
+        Input:
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+    
+        Return:
+            array of mx(degree + 1)
+    """
+    def get_normal_pars(self, x, y, m, seed = 1977):
+        
+        res = self.fit(x, y, method = "diff_evol")
+
+        assert res["success"], "Fit did not succeed."
+
+        # optimal parameters
+        pars = res["pars"]
+
+        # covariance matrix of parameters
+        cov = self.cov(x, y, pars)
+
+        # init random generator
+        rng = np.random.default_rng(seed)
+
+        return rng.multivariate_normal(pars, cov, size = m)
+
+    """
+        Generate m parameters via non-parametric bootstrapping with a minimal constraint 
+        that both groups should be present in the sampled data:
+
+        boostrapped sample = (xb, yb)  by sampling with replacement pairs (x_i, y_i) 
+                                    with condition that yb can not be just 0 or just 1 
+
+        Input:
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+    
+        Return:
+            array of mx(degree + 1)
+    """
+    def get_nonparam_boots_pars(self, x, y, m, seed = 1977):
+        
+        # fitting original data
+        res_fit = self.fit(x, y, method = "diff_evol")
+        assert res_fit["success"]
+
+        pars0 = res_fit["pars"]
+
+        n = len(x)
+
+        rng = np.random.default_rng(seed)
+
+        # generate parameters
+        lst = [pars0]
+        while True:
+
+            # create set indices for sampling with replacement + constraint
+            idx = rng.choice(n, n)
+            if np.sum(y[idx]) in [0, n]: continue
+
+            # do fitting 
+            res_fit = self.fit(x[idx], y[idx], pars0, method="local")
+            if not res_fit["success"]: continue
+            
+            # store pars
+            lst.append(res_fit["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)
+
+
+    """
+        Generate m parameters via non-parametric stratified bootstrapping:
+
+            boostrapped sample = (xb, yb)  
+                
+                xb = (sampled with replacement from x0, sampled with replacement from x1)
+                yb = (0...0, 1...1)
+
+        Input: 
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+        
+        Return:
+            array of mx(degree + 1)
+    """
+    def get_nonparam_stratified_boots_pars(self, x, y, m, seed = 1977):
+
+        # pars of original data
+        res_fit = self.fit(x, y, method = "diff_evol")
+        assert res_fit["success"]
+
+        pars0 = res_fit["pars"]
+        
+        # statistics about groups
+        xs = [x[y == i] for i in range(2)]
+        ns = [len(e) for e in xs]
+
+        # common vector states
+        yb = np.concatenate([np.full(ns[i], i) for i in range(2)])
+
+        rng = np.random.default_rng(seed)
+
+        # generate parameters
+        lst = [pars0]
+        for _ in range(m):
+            # stratified sampling with replacement
+            xb = np.concatenate([rng.choice(xs[i], ns[i]) for i in range(2)])
+
+            # do fitting 
+            res_fit = self.fit(xb, yb, pars0, method="local")
+            if not res_fit["success"]: continue
+            
+            # store pars
+            lst.append(res_fit["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)
+
+    """
+        Generate m parameters via parametric bootstrapping:
+
+            boostrapped sample = (x, yb) yb ~ B(model(x, fitted pars))  
+        
+        where B is Bernoulli distribution
+
+        Input: 
+            x : array of n floats
+            y : array of n int in {0,1}
+            m: integer, number of samples
+            seed : int, seed for the random generator
+        
+        Return:
+            array of mx(degree + 1)
+
+        Ref:
+
+        * https://www.scirp.org/journal/paperinformation?paperid=70962
+        * https://en.wikipedia.org/wiki/Bernoulli_distribution
+    """
+    def get_parametric_boots_pars(self, x, y, m, seed = 1977):
+
+        # first discuss original dataset
+        res_fit = self.fit(x, y, method = "diff_evol")
+        assert res_fit["success"]
+
+        pars0 = res_fit["pars"]
+
+        p = self.model(x, pars0)
+    
+        n = len(x)
+        
+        rng = np.random.default_rng(seed)
+
+        # generate parameters
+        lst = [pars0]
+        while True:
+
+            # Generate new binary outcomes from Bernoulli(p_i)
+            y_sim = rng.binomial(n = 1, p = p)
+            if np.sum(y_sim) in [0, n]: continue
+            
+            # do fitting 
+            res_fit = self.fit(x, y_sim, pars0, method="local")
+            if not res_fit["success"]: continue
+            
+            # store pars
+            lst.append(res_fit["pars"])
+            if len(lst) == m: break
+
+        return np.array(lst)

python/logistic/mono_cubic1.py

@@ -0,0 +1,120 @@
+# Function supporting monotonic cubic polynomials: first parametrization
+# We define the polynomial:
+#       poly(x) = sum_i beta_i * x^i
+# where the coefficients beta_i are parameterized by the vector 'pars'.
+
+import numpy as np
+
+def forward_map(pars, small=1e-8):
+    """
+    Map parameter vector pars to coefficients beta.
+
+    beta = [p0, (p2^2 + p1^2) / a, p2, p3^2],
+    where a = 3*p3^2 + small.
+    """
+    pars = np.asarray(pars, dtype=float)
+    p0, p1, p2, p3 = pars
+    a = 3.0 * p3**2 + small
+    beta = np.array([
+        p0,
+        (p2**2 + p1**2) / a,
+        p2,
+        p3**2
+    ], dtype=float)
+    return beta
+
+def backward_map(beta, small=1e-8):
+    """
+    Map coefficients beta back to parameter vector pars.
+    Nonnegative roots are used for p1 and p3.
+    """
+    beta = np.asarray(beta, dtype=float)
+    b0, b1, b2, b3 = beta
+    p0 = b0
+    p3 = np.sqrt(b3)
+    a = 3.0 * b3 + small
+    p2 = b2
+    p1_sq = b1 * a - p2**2
+    if p1_sq < 0:
+        raise ValueError(f"Negative square root encountered for p1^2 = {p1_sq}")
+    p1 = np.sqrt(p1_sq)
+    return np.array([p0, p1, p2, p3], dtype=float)
+
+def forward_map_jacobian(pars, small=1e-8):
+    """
+    Compute Jacobian of forward_map at given pars.
+    Returns J: d(beta)/d(pars), shape (4, 4).
+    """
+    p0, p1, p2, p3 = np.asarray(pars, dtype=float)
+    a = 3.0 * p3**2 + small
+    J = np.zeros((4, 4), dtype=float)
+
+    # beta0 row
+    J[0, 0] = 1.0
+
+    # beta1 row
+    J[1, 1] = 2.0 * p1 / a
+    J[1, 2] = 2.0 * p2 / a
+    J[1, 3] = -6.0 * p3 * (p2**2 + p1**2) / a**2
+
+    # beta2 row
+    J[2, 2] = 1.0
+
+    # beta3 row
+    J[3, 3] = 2.0 * p3
+
+    return J
+
+def forward_map_hessian(pars, small=1e-8):
+    """
+    Compute Hessians of forward_map at given pars.
+
+    Returns H: shape (4, 4, 4),
+    where H[i] is the 4x4 Hessian of beta[i] w.r.t. pars.
+    """
+    p0, p1, p2, p3 = np.asarray(pars, dtype=float)
+    a = 3.0 * p3**2 + small
+    N = p1**2 + p2**2
+    H = np.zeros((4, 4, 4), dtype=float)
+
+    # beta0: all zeros
+
+    # beta1 Hessian
+    H1 = np.zeros((4, 4), dtype=float)
+    H1[1, 1] = 2.0 / a
+    H1[2, 2] = 2.0 / a
+    H1[1, 3] = -12.0 * p1 * p3 / a**2
+    H1[3, 1] = H1[1, 3]
+    H1[2, 3] = -12.0 * p2 * p3 / a**2
+    H1[3, 2] = H1[2, 3]
+    H1[3, 3] = 6.0 * N * (12.0 * p3**2 - a) / a**3
+    H[1] = H1
+
+    # beta2: all zeros
+
+    # beta3 Hessian
+    H3 = np.zeros((4, 4), dtype=float)
+    H3[3, 3] = 2.0
+    H[3] = H3
+
+    return H
+
+
+# -------------------------
+# Round-trip test
+# -------------------------
+if __name__ == "__main__":
+    pars_original = np.array([1.0, 2.0, 3.0, 4.0])
+    beta = forward_map(pars_original)
+    pars_recovered = backward_map(beta)
+
+    print("Original pars: ", pars_original)
+    print("Beta:          ", beta)
+    print("Recovered pars:", pars_recovered)
+    print("Difference:    ", pars_recovered - pars_original)
+
+    print("\nJacobian at pars:")
+    print(forward_map_jacobian(pars_original))
+
+    print("\nHessian for beta1:")
+    print(forward_map_hessian(pars_original)[1])

python/logistic/mono_cubic2.py

@@ -0,0 +1,196 @@
+# Function supporting monotonic cubic polynomials: first parametrization
+# We define the polynomial:
+#       poly(x) = sum_i beta_i * x^i
+# where the coefficients beta_i are parameterized by the vector 'pars'.
+
+import numpy as np
+
+# Forward map: Parameters to polynomial coefficients
+def forward_map(pars):
+    """
+    Converts parameters (pars = [C, epsilon, k1, k2]) into polynomial coefficients (beta = [d, c, b, a]).
+    
+    Polynomial Definition:
+    - d = C: Constant term of the polynomial.
+    - c = k2^2 + epsilon^2: Coefficient of the linear term.
+    - b = k1 * k2: Coefficient of the quadratic term.
+    - a = k1**2 / 3: Coefficient of the cubic term.
+
+    Parameters:
+    - pars: NumPy array of parameters [C, epsilon, k1, k2].
+
+    Returns:
+    - beta: NumPy array of polynomial coefficients [d, c, b, a].
+    """
+    C, epsilon, k1, k2 = pars  # Unpack the parameter vector
+    
+    # Compute coefficients
+    d = C
+    c = k2**2 + epsilon**2
+    b = k1 * k2
+    a = k1**2 / 3
+    
+    # Return coefficients as a NumPy array
+    beta = np.array([d, c, b, a])
+    return beta
+
+
+# Jacobian of the forward map: First-order derivatives
+def forward_map_jacobian(pars):
+    """
+    Computes the Jacobian matrix of the forward map analytically.
+    
+    Parameters:
+    - pars: NumPy array of parameters [C, epsilon, k1, k2].
+    
+    Returns:
+    - J: NumPy 4x4 Jacobian matrix, where J[i, j] = d(beta[i])/d(pars[j]).
+    """
+    C, epsilon, k1, k2 = pars  # Unpack parameters
+    
+    # Initialize Jacobian matrix
+    J = np.zeros((4, 4))  # 4x4 matrix
+    
+    # Partial derivatives for d = C
+    J[0, 0] = 1  # d(d)/dC
+    J[0, 1] = 0  # d(d)/d(epsilon)
+    J[0, 2] = 0  # d(d)/d(k1)
+    J[0, 3] = 0  # d(d)/d(k2)
+
+    # Partial derivatives for c = k2^2 + epsilon^2
+    J[1, 0] = 0  # d(c)/dC
+    J[1, 1] = 2 * epsilon  # d(c)/d(epsilon)
+    J[1, 2] = 0  # d(c)/d(k1)
+    J[1, 3] = 2 * k2  # d(c)/d(k2)
+
+    # Partial derivatives for b = k1 * k2
+    J[2, 0] = 0  # d(b)/dC
+    J[2, 1] = 0  # d(b)/d(epsilon)
+    J[2, 2] = k2  # d(b)/d(k1)
+    J[2, 3] = k1  # d(b)/d(k2)
+
+    # Partial derivatives for a = k1^2 / 3
+    J[3, 0] = 0  # d(a)/dC
+    J[3, 1] = 0  # d(a)/d(epsilon)
+    J[3, 2] = 2 * k1 / 3  # d(a)/d(k1)
+    J[3, 3] = 0  # d(a)/d(k2)
+    
+    return J
+
+
+# Hessian of the forward map: Second-order derivatives
+def forward_map_hessian(pars):
+    """
+    Computes the Hessian tensor of the forward map analytically.
+    
+    Parameters:
+    - pars: NumPy array of parameters [C, epsilon, k1, k2].
+    
+    Returns:
+    - H: NumPy 4x4x4 Hessian tensor, where H[i, j, k] = d^2(beta[i])/d(pars[j])d(pars[k]).
+    """
+    C, epsilon, k1, k2 = pars  # Unpack parameters
+    
+    # Initialize Hessian tensor (4 x 4 x 4)
+    H = np.zeros((4, 4, 4))
+    
+    # Hessian for d = C: All second derivatives are zero
+    # Already H[0, :, :] is initialized to zero
+    
+    # Hessian for c = k2^2 + epsilon^2
+    H[1, 1, 1] = 2  # d^2(c)/d(epsilon^2)
+    H[1, 3, 3] = 2  # d^2(c)/d(k2^2)
+    
+    # Hessian for b = k1 * k2: All second derivatives are zero
+    # Already H[2, :, :] is initialized to zero
+    
+    # Hessian for a = k1^2 / 3
+    H[3, 2, 2] = 2 / 3  # d^2(a)/d(k1^2)
+    
+    return H
+
+
+# Backward map: Polynomial coefficients to parameters
+def backward_map(beta, only_one = True):
+    """
+    Computes the parameters (pars = [C, epsilon, k1, k2]) from the polynomial coefficients (beta = [d, c, b, a]).
+    
+    Polynomial Definition:
+    - d = C: Constant term of the polynomial.
+    - c = k2^2 + epsilon^2: Used to recover k2 and epsilon.
+    - b = k1 * k2: Used to recover k1 and k2.
+    - a = k1^2 / 3: Used to recover k1.
+
+    Parameters:
+    - beta: NumPy array of polynomial coefficients [d, c, b, a].
+
+    Returns:
+    - List of possible parameter sets [(C, epsilon, k1, k2)].
+    """
+    d, c, b, a = beta  # Unpack the coefficients
+    
+    # Recover k1 from a (two possible values due to ± sqrt)
+    if a < 0:
+        raise ValueError("Coefficient 'a' must be non-negative for monotonic polynomials.")
+    
+    k1_options = np.unique([np.sqrt(3 * a), -np.sqrt(3 * a)])  # Two possible k1 values
+    
+    possible_parameters = []
+    
+    # For each possible k1, compute k2 and epsilon
+    for k1 in k1_options:
+        k2 = None
+        if k1 != 0:  # Ensure k1 is non-zero (avoids division by zero)
+            k2 = b / k1  # Compute k2 from b and k1
+        elif b == 0:
+            k2 = 0
+
+        if k2 is None: continue    
+       
+        # Check if c >= k2^2 for valid epsilon computation
+        if c >= k2**2:
+            epsilon_options = np.unique([np.sqrt(c - k2**2), -np.sqrt(c - k2**2)])  # Two possible epsilon values
+            
+            for epsilon in epsilon_options:
+                # Constant term d maps directly to C
+                C = d
+                sol = np.array([C, epsilon, k1, k2])
+
+                if only_one: return sol
+                possible_parameters.append(sol)
+
+
+    return possible_parameters
+
+
+# -------------------------
+# Round-trip test
+# -------------------------
+if __name__ == "__main__":
+    pars_original = np.array([1.0, 2.0, 3.0, 4.0])
+    beta = forward_map(pars_original)
+    pars_recovered = backward_map(beta)
+
+    print("Original pars: ", pars_original)
+    print("Beta:          ", beta)
+    print("Recovered pars:", pars_recovered)
+    print("Difference:    ", pars_recovered - pars_original)
+
+    print("\nJacobian at pars:")
+    print(forward_map_jacobian(pars_original))
+
+    print("\nHessian for beta1:")
+    print(forward_map_hessian(pars_original)[1])
+
+    print("\nLinear func:")
+    lin_fun_beta = [1,0.2,0,0]
+    only_one = False
+
+    lin_fun_pars = backward_map(lin_fun_beta, only_one=only_one)
+    lin_fun_beta_recover = lin_fun_pars if only_one else np.unique([forward_map(pars) for pars in lin_fun_pars], axis=0)
+    
+    print(f"  {only_one = }")
+    print("  lin_fun_beta:", lin_fun_beta)
+    print("  backwards:", lin_fun_pars)
+    print("  forwards:", lin_fun_beta_recover)
+

python/logistic/zahra_plus_comments.py

@@ -0,0 +1,141 @@
+"""
+    Working on Bayesian formula model + penalty  == MAP approach with prior = penalty
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import optimize, stats
+
+#here is data
+assert np.all(X > 0), "All X must be > 0"
+n1, n0 = int(y.sum()), int((1-y).sum())
+
+m_emp = n1 / (n1 + n0)  # approx Prob(AE)
+
+print(f"AE={n1}, NC={n0}, empirical p_AE={m_emp:.6f}")
+
+# Model: AE ~ BetaPrime(a,b,scale), NC ~ LogNormal(mu, sigma)
+# p_ae = Prob(AE)
+# θ = [p_ae, a, b, sc, mu_nc, sig_nc]
+def logpdf_betaprime(x, a, b, scale):
+    return stats.betaprime.logpdf(x, a=a, b=b, scale=scale)
+
+def logpdf_lognorm(x, mu, sigma):
+    return stats.lognorm.logpdf(x, s=sigma, scale=np.exp(mu))
+
+def neg_conditional_ll(theta, X, y, eps=1e-12):
+    
+    p_ae, a, b, sc, mu_nc, sig_nc = theta
+
+    logf1 = logpdf_betaprime(X, a, b, sc)     # AE
+    logf0 = logpdf_lognorm(X, mu_nc, sig_nc)  # NC
+
+    # computing p(AE|x) = 1/(1 + exp(-logit))
+    logit = np.log(p_ae) - np.log(1.0 - p_ae) + (logf1 - logf0)
+    p = 1.0 / (1.0 + np.exp(-np.clip(logit, -50, 50)))
+
+    # computing general nllf = -llf
+    #  llf = sum_i log(p(y_i|x_i)) 
+    #      = sum_i y_i log( p(AE|x) + (1- y_i) log(1 - p(AE|x));  p(NC|x) = 1 -p(AE|x)
+
+    nllf = -np.sum(y*np.log(p+eps) + (1-y)*np.log(1-p+eps))
+    return nllf
+
+"""
+    MAP regularization
+        mean = m_emp, concentration = TAU
+    using beta distribution B(alpha, beta). The parameters are set as
+        alpha = m_emp * TAU, beta = (1-m_emp) * TAU  
+    this yields E[X] = alpha/(alpha+beta) = m_emp
+        Var[X] = m_emp(1- m_emp)/(tau +1)  
+    this could be variance between institutions collected by Katja, 
+    my estimate is 
+        var = 5%^2 
+    This yields tau about 30. 
+"""
+TAU = 100.0            # ↑ increase to pull p_AE closer to empirical prior
+alpha = max(m_emp * TAU, 1e-6)
+beta  = max((1 - m_emp) * TAU, 1e-6)
+
+
+def neg_log_prior_p(p, eps=1e-12):
+    # -log Beta(p | alpha, beta) up to a constant
+    return -( (alpha - 1)*np.log(p + eps) + (beta - 1)*np.log(1 - p + eps) )
+
+
+def neg_posterior(theta, X, y):
+    nll = neg_conditional_ll(theta, X, y)
+    return nll + neg_log_prior_p(theta[0])   # add prior penalty on p_AE only
+
+
+# Bounds (optionally enforce heavy AE tail with b ≤ 1)
+# =========================
+HEAVY_TAIL = True     # set False if you don't want to force the right asymptote
+b_upper = 1.0 if HEAVY_TAIL else 50.0
+bounds = [
+    (1e-3, 1-1e-3),                           # p_ae (free, but regularized by Beta prior)
+    (0.20, 50.0),                             # a (AE BetaPrime)
+    (0.20, b_upper),                          # b (AE BetaPrime)  <-- heavy tail if ≤ 1
+    (0.01, 10.0),                             # scale (AE BetaPrime)
+    (np.log(X).min()-2.0, np.log(X).max()+2.0),  # mu_nc (LogNormal)
+    (0.05, 2.0),                              # sigma_nc (LogNormal)
+]
+# Initialization
+p0 = np.clip(m_emp, bounds[0][0], bounds[0][1])
+X0 = X[y==0]
+mu0 = float(np.mean(np.log(X0)))
+sig0 = float(np.std(np.log(X0), ddof=0))
+mu0 = np.clip(mu0, bounds[4][0], bounds[4][1])
+sig0 = np.clip(sig0, bounds[5][0], bounds[5][1])
+X1 = X[y==1]
+m1 = float(np.mean(X1))
+a0, b0 = 2.5, min(0.8, b_upper)  # start with heavy-tail-ish b if allowed
+sc0 = np.clip(m1 * (b0 - 1 + 1e-6) / max(a0, 1e-6), bounds[3][0], bounds[3][1])
+theta0 = np.array([p0, a0, b0, sc0, mu0, sig0], dtype=float)
+# FREE prior (for comparison)
+# =========================
+res_free = optimize.minimize(
+    fun=neg_conditional_ll,
+    x0=theta0,
+    args=(X, y),
+    method="L-BFGS-B",
+    bounds=bounds,
+    options=dict(maxiter=4000, ftol=1e-12)
+)
+theta_free = res_free.x
+print("\n[FREE prior] p_AE =", float(theta_free[0]), "   CLL =", -res_free.fun)
+# MAP prior (regularized toward empirical)
+# =========================
+res_map = optimize.minimize(
+    fun=neg_posterior,
+    x0=theta0,
+    args=(X, y),
+    method="L-BFGS-B",
+    bounds=bounds,
+    options=dict(maxiter=4000, ftol=1e-12)
+)
+theta_map = res_map.x
+print("[MAP prior]  p_AE =", float(theta_map[0]), "   CLL(post) =", -res_map.fun)
+# MAP parameters
+p_hat, a_hat, b_hat, sc_hat, mu_hat, sig_hat = theta_map
+print("\nFitted (MAP) parameters:")
+print(f"  p_AE = {p_hat:.6f}  (empirical {m_emp:.6f}, TAU={TAU})")
+print(f"  AE BetaPrime: a={a_hat:.4f}, b={b_hat:.4f}, scale={sc_hat:.4f}")
+print(f"  NC LogNormal: mu={mu_hat:.4f}, sigma={sig_hat:.4f}")
+# Posterior & plot
+def predict_proba(x):
+    x = np.asarray(x, dtype=float)
+    logf1 = logpdf_betaprime(x, a_hat, b_hat, sc_hat)
+    logf0 = logpdf_lognorm(x, mu_hat, sig_hat)
+    logit = np.log(p_hat) - np.log(1.0 - p_hat) + (logf1 - logf0)
+    return 1.0 / (1.0 + np.exp(-np.clip(logit, -50, 50)))
+x_grid = np.linspace(max(1e-6, X.min()*0.6), max(X.max()*2.0, 8.0), 600)
+p_grid = predict_proba(x_grid)
+plt.figure(figsize=(8,5))
+plt.plot(x_grid, p_grid, 'k-', linewidth=2, label="P(AE | X) [MAP]")
+plt.scatter(X[y==1], np.ones(n1), marker='x', label="AE samples")
+plt.scatter(X[y==0], np.zeros(n0), marker='o', label="NC samples")
+plt.xlabel("SUV feature X"); plt.ylabel("Predicted P(AE | X)")
+plt.title("BetaPrime–LogNormal with MAP prior on p_AE")
+plt.ylim(-0.05, 1.05); plt.legend(); plt.grid(True)
+plt.show() (edited) 

refs/references.txt

@@ -19,3 +19,7 @@ https://python-bloggers.com/2024/03/estimating-logistic-regression-coefficients-
 https://www.mathworks.com/help/stats/mlecov.html 
 
 
+Reference For Bayesian
+http://ndl.ethernet.edu.et/bitstream/123456789/37609/1/William%20M.%20Bolstad_2017.pdf
+
+