|
|
@@ -6,8 +6,6 @@ from scipy import optimize
|
|
|
from scipy.special import betaln, gammaln
|
|
|
import scipy.io as io
|
|
|
|
|
|
-# ---------- Data ----------
|
|
|
-data_path = "../data/"
|
|
|
|
|
|
# === Load data ===
|
|
|
data_path = "../data/"
|
|
|
@@ -25,7 +23,7 @@ X = np.concatenate([X_AE, X_NC], axis=0)
|
|
|
n = len(y); n1 = int(y.sum()); p_emp = n1 / n
|
|
|
rng = np.random.default_rng(12345)
|
|
|
|
|
|
-# ---------- Helpers ----------
|
|
|
+# Helpers
|
|
|
def logistic(z):
|
|
|
z = np.clip(z, -60, 60)
|
|
|
return 1.0/(1.0+np.exp(-z))
|
|
|
@@ -37,7 +35,7 @@ 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) ---
|
|
|
+# 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):
|
|
|
@@ -55,7 +53,7 @@ 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 ----------
|
|
|
+# Global monotonicity cap for theta
|
|
|
def theta_max(a, b, k, s, eps=1e-12):
|
|
|
A = a - k
|
|
|
if A <= 0:
|
|
|
@@ -92,7 +90,7 @@ def unpack_phi_mono(phi):
|
|
|
|
|
|
return p, a, b, s, k, th
|
|
|
|
|
|
-# ---------- Priors ----------
|
|
|
+# 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)
|
|
|
@@ -105,7 +103,7 @@ def nlog_lognormal(x, mu, sigma, eps=1e-12):
|
|
|
lx = np.log(x)
|
|
|
return 0.5 * ((lx - mu)/sigma)**2 + lx
|
|
|
|
|
|
-# ---------- Objective ----------
|
|
|
+# 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
|
|
|
@@ -118,7 +116,7 @@ def neg_post_phi_mono_WITH_CONST_REG(phi, X, y):
|
|
|
# Prior on p ~ Beta(alpha, beta)
|
|
|
npr_p = -((alpha-1)*np.log(p + eps) + (beta-1)*np.log(1 - p + eps))
|
|
|
|
|
|
- # --- Regularization (weak priors) ---
|
|
|
+ # Regularization (weak priors)
|
|
|
# AE median m1: use AE median if present; otherwise overall median.
|
|
|
if (y == 1).any():
|
|
|
m1 = np.median(X[y == 1])
|
|
|
@@ -150,7 +148,7 @@ def neg_post_phi_mono_WITH_CONST_REG(phi, X, y):
|
|
|
|
|
|
return nll + npr_p + reg + npr_r
|
|
|
|
|
|
-# ---------- Initialization ----------
|
|
|
+# 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
|
|
|
@@ -200,7 +198,7 @@ def init_phi(X, y):
|
|
|
], float)
|
|
|
return raw
|
|
|
|
|
|
-# ---------- Fitting ----------
|
|
|
+# 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)
|
|
|
@@ -220,7 +218,7 @@ def fit_hard_mono_WITH_CONST_REG(X, y, phi_start=None, maxtries=6, jitter=0.3, r
|
|
|
phi = phi + rng.normal(0, jitter, size=phi.shape)
|
|
|
raise RuntimeError(f"Fit failed. Last status: {getattr(last_err, 'message', 'n/a')}")
|
|
|
|
|
|
-# ---------- Convenience ----------
|
|
|
+# 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)
|
|
|
@@ -266,7 +264,7 @@ def plot_s_shape(theta, X, y, rng=None, ax=None, label='P(AE | x)'):
|
|
|
ax.legend(loc='lower right', frameon=False)
|
|
|
return ax
|
|
|
|
|
|
-# ---------- Run fit ----------
|
|
|
+# 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)
|
|
|
@@ -281,7 +279,7 @@ plt.show()
|
|
|
|
|
|
|
|
|
|
|
|
-# Delta-method
|
|
|
+# 1- Delta-method
|
|
|
import numdifftools as nd
|
|
|
|
|
|
# wrap scalar objective for numdifftools
|
