Cleaning

notes/20250606_minutes.txt

@@ -0,0 +1,71 @@
+Date: 8.june 2025
+
+Meeting with Zahra
+
+Dates:
+  30.6. - Zahra has presentation
+  5.7.  - Zahra visits Iran for 2-3 weeks
+
+Asking about accuracy measures:
+  * for fit we have goodness of fit measures:
+    - log likelihood
+    - AIC and BIC
+  * these models used for classification: [1]
+  -  pred Y = [ Y : f(x|Y) > 0.5 === decision function(x) > 0
+              [ !Y : otherwise
+  - in logit models: [2]  Speech and Language Processing, Ch 5
+      decision_function(x) = b + w^T x == logit (x)
+      logit == log odds
+      logit(x) = 1/(1 + exp(-decision_function(x))
+  - ACC == accuracy score [3]
+      accuracy(y, \hat{y}) = 1/n \sum_{i=0}^{n-1} delta_{\hat{y}_i, y_i}
+
+  - AUC == area under curve ???? [4,5]
+      From [6]:
+
+        fpr, tpr, thresholds = ROC = function(y_true, y_score)
+
+        We are systematically doing
+
+          pred Y[i] = [ 1 : y_score[i] >= threshold
+                     [ 0 : otherwise
+          all i
+
+        and comparing with
+
+          y_true[i] all i
+        for all meaningful thresholds == y_score
+
+      Target scores, can either be probability estimates of the positive class,
+      confidence values, or non-thresholded measure of decisions (as returned by
+      “decision_function” on some classifiers). For decision_function scores,
+      values greater than or equal to zero should indicate the positive class.
+      In our case y_score = prob estimate(for Y=1) as given by the model
+
+
+Ref:
+   - [1] https://realpython.com/logistic-regression-python/
+   - [2] https://web.stanford.edu/~jurafsky/slp3/5.pdf
+   - https://en.wikipedia.org/wiki/Logistic_regression
+
+   - [3] https://scikit-learn.org/stable/modules/model_evaluation.html#accuracy-score
+   - [4] https://scikit-learn.org/stable/modules/model_evaluation.html#roc-metrics
+   - [5] https://en.wikipedia.org/wiki/Receiver_operating_characteristic
+   - [6] https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_curve.html
+
+   - 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
+
+Plan:
+
+  * Make notes overleaf file (ZA)
+    - notes will the basis for the paper and presentations
+  * Transfer as much text and results to notes (ZA)
+  * Test calc. with bayesian model (MH)
+  * Chat about progress at the picnic
+
+Ideal:
+  * Have some preliminary results for bayesian model until Friday
+
+

notes/PipeLine.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "id": "4332243f",
    "metadata": {},
    "outputs": [
@@ -35,8 +35,7 @@
     }
    ],
    "source": [
-    "# condprob_modeling.py\n",
-    "\n",
+    "#PipleLine for Conditional Probability Modeling\n",
     "def define_model():\n",
     "    print(\"[Task 1] Define a model for conditional probability P(X|Y)\")\n",
     "    print(\"        - Specify assumptions\")\n",

notes/Untitled-1.ipynb

@@ -0,0 +1,24 @@
+{
+ "cells": [],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "base",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

python/# BetaPrime vs Gamma, hard-mono fit WITH.py

@@ -0,0 +1,187 @@
+# 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 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/Untitled-1.ipynb

@@ -0,0 +1,208 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eae6df92",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Constrained Bayesian fit (Gamma for NC, Beta-Prime for AE) — no regularization, no r-prior\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from scipy import optimize\n",
+    "from scipy.special import betaln, gammaln\n",
+    "import scipy.io as io\n",
+    "\n",
+    "data_path = \"../data/\"\n",
+    "suv   = io.loadmat(data_path + \"suv_percentilesSLOthenUWM.mat\")['lung_SUVperc_COMBINED'][0:58, :, :]\n",
+    "flags = io.loadmat(data_path + \"flags_combined.mat\")['flags'][0:58, 3]   # 0=NC, 1=AE\n",
+    "\n",
+    "# Feature X = max SUV_94 per subject; label y = flags\n",
+    "X = np.nanmax(suv[:, :, 94], axis=1).astype(float).ravel()\n",
+    "y = np.asarray(flags, int).ravel()\n",
+    "\n",
+    "# Guard for logs\n",
+    "X = np.clip(X, 1e-12, None)\n",
+    "p_emp = float(y.mean())\n",
+    "\n",
+    "# Small helpers\n",
+    "\n",
+    "def logistic(z):\n",
+    "    z = np.clip(z, -60, 60)\n",
+    "    return 1.0 / (1.0 + np.exp(-z))\n",
+    "\n",
+    "def sigmoid(t):\n",
+    "    return 1.0 / (1.0 + np.exp(-t))\n",
+    "\n",
+    "def softplus(t):\n",
+    "    t = np.asarray(t, float)\n",
+    "    return np.log1p(np.exp(-np.abs(t))) + np.maximum(t, 0.0)\n",
+    "\n",
+    "# dE(x) pieces for log-odds\n",
+    "def dE_ess(x, a, b, s, k, th):\n",
+    "    x = np.asarray(x, float)\n",
+    "    return (a - k) * np.log(x) - (a + b) * np.log1p(x / s) + x / th\n",
+    "\n",
+    "def dE_const(a, b, s, k, th):\n",
+    "    return -(a * np.log(s)) - betaln(a, b) + k * np.log(th) + gammaln(k)\n",
+    "\n",
+    "def dE_full(x, a, b, s, k, th):\n",
+    "    return dE_ess(x, a, b, s, k, th) + dE_const(a, b, s, k, th)\n",
+    "\n",
+    "# Monotonicity cap for theta\n",
+    "def theta_max(a, b, k, s, eps=1e-12):\n",
+    "    A = a - k\n",
+    "    if A <= 0:\n",
+    "        return np.inf\n",
+    "    r = np.sqrt(a + b) - np.sqrt(max(A, eps))\n",
+    "    if r <= 1e-12:\n",
+    "        return np.inf\n",
+    "    return s / (r * r)\n",
+    "\n",
+    "# φ = [p_raw, b_raw, s_raw, k_raw, d_raw, u_raw]  (unconstrained)\n",
+    "def unpack_phi_mono(phi):\n",
+    "    p_raw, b_raw, s_raw, k_raw, d_raw, u_raw = phi\n",
+    "    p = sigmoid(p_raw)                       # (0,1)\n",
+    "    b = softplus(b_raw) + 1e-6               # >0\n",
+    "    s = softplus(s_raw) + 1e-6               # >0\n",
+    "    k = softplus(k_raw) + 1e-6               # >0\n",
+    "    delta = softplus(d_raw) + 1e-6           # >0\n",
+    "    a = k + delta                            # enforce a > k\n",
+    "    th_cap = theta_max(a, b, k, s)           # theta cap from monotonicity\n",
+    "    th = th_cap * sigmoid(u_raw)             # 0 < theta <= th_cap\n",
+    "    return p, a, b, s, k, th\n",
+    "\n",
+    "#  Prior on p \n",
+    "TAU = 25.0                                   # shrink toward empirical AE rate\n",
+    "alpha = max(TAU * p_emp, 1e-6)\n",
+    "beta  = max(TAU * (1.0 - p_emp), 1e-6)\n",
+    "\n",
+    "prior_r = None\n",
+    "#prior_r = (1.01, 1.01)\n",
+    "#prior_r = (1.05, 1.05)\n",
+    "#prior_r = (3, 3)\n",
+    "#prior_r = (1.2, 1.2)\n",
+    "# Objective: negative log-posterior (likelihood + Beta prior on p)\n",
+    "def neg_post_phi_mono(phi, X, y):\n",
+    "    p, a, b, s, k, th = unpack_phi_mono(phi)\n",
+    "    eps = 1e-12\n",
+    "\n",
+    "    L  = (np.log(p) - np.log(1 - p)) + dE_full(X, a, b, s, k, th)\n",
+    "    px = logistic(L)\n",
+    "    nll = -np.sum(y * np.log(px + eps) + (1 - y) * np.log(1 - px + eps))\n",
+    "\n",
+    "    # Beta(alpha, beta) prior on p → negative log-prior\n",
+    "    npr_p = -((alpha - 1) * np.log(p + eps) + (beta - 1) * np.log(1 - p + eps))\n",
+    "\n",
+    "    if prior_r is None: return nll + npr_p\n",
+    "        \n",
+    "    # Prior on r = theta / theta_max (softly avoid boundaries)\n",
+    "    thcap = theta_max(a, b, k, s)\n",
+    "    if np.isfinite(thcap) and thcap > 0:\n",
+    "        r = np.clip(th / thcap, 1e-9, 1 - 1e-9)\n",
+    "        # Negative log Beta prior:  -[(α-1)log r + (β-1)log(1-r)]  (const dropped)\n",
+    "        npr_r = -((prior_r[0] -1) * np.log(r) + (prior_r[1]-1)* np.log(1.0 - r))\n",
+    "    else:\n",
+    "        npr_r = 0.0\n",
+    "\n",
+    "    return nll + npr_p + npr_r\n",
+    "\n",
+    "# Initialization (stable, simple)\n",
+    "def init_phi(X, y):\n",
+    "    # Gamma(k, theta) MoM for NC group (only need k0 as a safe size proxy)\n",
+    "    X0 = X[y == 0]\n",
+    "    m0 = X0.mean() if X0.size else X.mean()\n",
+    "    v0 = X0.var()  if X0.size else X.var()\n",
+    "    k0 = 2.0 if v0 <= 0 else max((m0**2)/(v0 + 1e-9), 1.5)\n",
+    "\n",
+    "    # AE median to seed s0\n",
+    "    X1 = X[y == 1]\n",
+    "    m1 = np.median(X1) if X1.size else np.median(X)\n",
+    "\n",
+    "    p0 = np.clip(float(y.mean()), 1e-3, 1 - 1e-3)\n",
+    "    b0, s0 = 1.5, max(m1, 0.5)\n",
+    "\n",
+    "    return np.array([\n",
+    "        np.log(p0 / (1 - p0)),             # p_raw\n",
+    "        np.log(np.expm1(b0) + 1e-9),       # b_raw\n",
+    "        np.log(np.expm1(s0) + 1e-9),       # s_raw\n",
+    "        np.log(np.expm1(k0) + 1e-9),       # k_raw\n",
+    "        np.log(np.expm1(1.0) + 1e-9),      # d_raw  (delta)\n",
+    "        -0.2                               # u_raw  (keeps theta a bit below cap initially)\n",
+    "    ], float)\n",
+    "\n",
+    "#  Fit wrapper (one retry)\n",
+    "def fit_bayes_mono(X, y, phi_start=None, rng=None):\n",
+    "    if rng is None:\n",
+    "        rng = np.random.default_rng(0)\n",
+    "    if phi_start is None:\n",
+    "        phi_start = init_phi(X, y)\n",
+    "\n",
+    "    obj = lambda phi: neg_post_phi_mono(phi, X, y)\n",
+    "\n",
+    "    res = optimize.minimize(\n",
+    "        obj, phi_start, method=\"L-BFGS-B\",\n",
+    "        options={\"maxiter\": 6000, \"ftol\": 1e-9}\n",
+    "    )\n",
+    "    if not (res.success and np.isfinite(res.fun)):\n",
+    "        phi_try = phi_start + rng.normal(0, 0.2, size=phi_start.shape)\n",
+    "        res = optimize.minimize(\n",
+    "            obj, phi_try, method=\"L-BFGS-B\",\n",
+    "            options={\"maxiter\": 6000, \"ftol\": 1e-9}\n",
+    "        )\n",
+    "    return unpack_phi_mono(res.x), res\n",
+    "\n",
+    "#  prediction \n",
+    "def P_with(theta, x):\n",
+    "    p, a, b, s, k, th = theta\n",
+    "    L = (np.log(p) - np.log(1 - p)) + dE_full(x, a, b, s, k, th)\n",
+    "    return logistic(L)\n",
+    "\n",
+    "# Run fit + plot\n",
+    "theta_hat, res = fit_bayes_mono(X, y)\n",
+    "print(\"Optimization success:\", res.success, \" fval:\", float(res.fun))\n",
+    "\n",
+    "(p,a,b,s,k,th) = theta_hat\n",
+    "thcap = theta_max(a, b, k, s)\n",
+    "print(\"theta (p,a,b,s,k,theta):\", tuple(float(t) for t in theta_hat), \"ratio(th):\", th/thcap)\n",
+    "\n",
+    "# x-range (cap right end at 10 for readability)\n",
+    "x_lo = max(1e-6, float(X.min()) * 0.8)\n",
+    "x_hi = min(10.0, float(X.max()) * 1.2)\n",
+    "xg   = np.linspace(x_lo, x_hi, 600)\n",
+    "p_curve = P_with(theta_hat, xg)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(7.0, 4.6), dpi=140)\n",
+    "ax.plot(xg, p_curve, color=\"#000000\", lw=2.2, label=\"P(AE|x) (MAP)\")\n",
+    "\n",
+    "# overlay data with tiny vertical jitter so points don't overlap\n",
+    "rng_plot = np.random.default_rng(999)\n",
+    "jit = (rng_plot.random(len(y)) - 0.5) * 0.06\n",
+    "ax.scatter(X[y==0], (y + jit)[y==0], s=22, alpha=0.55, color=\"#2ca02c\", edgecolors='none', label='NC')\n",
+    "ax.scatter(X[y==1], (y + jit)[y==1], s=26, alpha=0.75, color=\"#ff7f0e\", edgecolors='none', label='AE')\n",
+    "\n",
+    "ax.set_ylim(-0.05, 1.05)\n",
+    "ax.set_xlabel('x')\n",
+    "ax.set_ylabel('P(AE | x)')\n",
+    "\n",
+    "if prior_r is None:\n",
+    "    ax.set_title('Constrained Bayesian fit (no regularization, prior on p)')\n",
+    "else:\n",
+    "    ax.set_title(f'Constrained Bayesian fit (no regularization, prior on p and prior r{prior_r})')\n",
+    "\n",
+    "ax.grid(alpha=0.3)\n",
+    "ax.legend(loc='lower right', frameon=False)\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}