Update 'python/bayesian/Bayesian_Zahra.py'

python/bayesian/Bayesian_Zahra.py

@@ -280,78 +280,65 @@ plt.show()
 
 
 # 1- 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
-
+# ===== Delta band + compact summaries (minimal) =====
+import numpy as np
+import matplotlib.pyplot as plt
+import numdifftools as nd
+from scipy.stats import norm
 
-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
+# 1) Covariance in raw-phi space at MAP
+phi_hat   = res.x.copy()
+f_obj     = lambda phi: neg_post_phi_mono_WITH_CONST_REG(np.asarray(phi, float), X, y)
+H         = nd.Hessian(f_obj, method='central')(phi_hat)
+Sigma_phi = np.linalg.pinv(0.5*(H + H.T))  # robust inverse
 
-# 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)
+# 2) Delta band on P(AE|x) via numdifftools.Gradient
+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%
+p_hat = np.empty_like(xg)
+p_lo  = np.empty_like(xg)
+p_hi  = np.empty_like(xg)
 
 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))
-
+    gx   = g_px(x)
+    ph   = gx(phi_hat)
+    grad = nd.Gradient(gx, method='central')(phi_hat)
+    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()
-
-
+    p_lo[i]  = np.clip(ph - z*se, 0.0, 1.0)
+    p_hi[i]  = np.clip(ph + z*se, 0.0, 1.0)
+
+# 3) Plot
+fig, ax = plt.subplots(figsize=(7.2, 4.4), dpi=140)
+ax.plot(xg, p_hat, lw=2.0, label='P(AE|x) @ MAP')
+ax.fill_between(xg, p_lo, p_hi, alpha=0.20, label='95% Delta band')
+rngp = np.random.default_rng(999); jit = (rngp.random(len(X)) - 0.5) * 0.06
+ax.scatter(X[y==0], (y+jit)[y==0], s=22, alpha=0.55, edgecolors='none', label='NC')
+ax.scatter(X[y==1], (y+jit)[y==1], s=26, alpha=0.75, edgecolors='none', label='AE')
+ax.set_ylim(-0.05, 1.05); ax.set_xlabel('x'); ax.set_ylabel('P(AE | x)')
+ax.grid(alpha=0.3); ax.legend(loc='lower right')
+plt.tight_layout(); plt.show()
+
+
+
+w = p_hi - p_lo
+mask = (xg >= float(X.min())) & (xg <= float(X.max()))
+print("\nBand width (95% pointwise): "
+      f"overall mean {w.mean():.3f}, max {w.max():.3f}; "
+      f"in-range mean {w[mask].mean():.3f}, max {w[mask].max():.3f}")
+
+
+try:
+    G = lambda phi: np.array(unpack_phi_mono(np.asarray(phi, float)), float)  # -> [p,a,b,s,k,theta]
+    J = nd.Jacobian(G)(phi_hat)
+    Sigma_theta = J @ Sigma_phi @ J.T
+    se = np.sqrt(np.maximum(np.diag(Sigma_theta), 0.0))
+    theta_hat_vec = G(phi_hat); names = ["p","a","b","s","k","theta"]
+    print("\nParameter 95% CIs (Delta/Wald):")
+    for nm, v, svi in zip(names, theta_hat_vec, se):
+        print(f"  {nm:>6s} : {v:.6g}  [ {v - z*svi:.6g}, {v + z*svi:.6g} ]")
+except Exception as e:
+    print("(Parameter CI step skipped:", e, ")")