Thesis - wip

This commit is contained in:
2026-07-15 20:58:17 +03:00
parent ce726d938f
commit a85bb2a551
13 changed files with 963 additions and 399 deletions
Binary file not shown.

After

Width:  |  Height:  |  Size: 76 KiB

@@ -0,0 +1,181 @@
"""
Comparison of density estimators on real model output scores.
For each class (negative / positive) plots four curves against the empirical CDF:
- ECDF — ground truth step function
- KDE (ISJ) — Gaussian-kernel KDE with Improved Sheather-Jones bandwidth
- Gaussian fit — MLE normal distribution
- Beta fit — sigmoid-transformed scores fitted with Beta MLE, CDF
mapped back to the original logit axis
Data: test split of a ROLL model trained on the glass0 dataset.
Run from the impl directory (KDEpy lives in that flake):
cd impl && nix develop --command python3 ../thesis/content/method/figures/density_fit_comparison.py
"""
import sys, pickle
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import norm, beta as sp_beta
from scipy.special import expit
from KDEpy import FFTKDE
from KDEpy.bw_selection import silvermans_rule, improved_sheather_jones
# Resolve paths relative to this script's location so the script runs from anywhere.
_HERE = Path(__file__).resolve().parent # .../thesis/content/method/figures
_THESIS = _HERE.parents[2] # .../thesis
_IMPL = _THESIS.parent / "impl" # .../impl
sys.path.insert(0, str(_IMPL)) # pickle needs 'import src.experiment' → impl/src/experiment.py
# ── Config ──────────────────────────────────────────────────────────────────
RESULT_PKL = _IMPL / "results/glass0/2026-07-08-22-25/roll/0/test-res.pkl"
OUT_DIR = _HERE
SPLIT = "test"
C_NEG = "#4477AA"
C_POS = "#CC6633"
plt.rcParams.update({
"text.usetex": False,
"mathtext.fontset": "cm",
"pdf.fonttype": 42,
})
# ── Load data ────────────────────────────────────────────────────────────────
with open(RESULT_PKL, "rb") as f:
ep = pickle.load(f)
res = ep.split_results[SPLIT]
yh = res.yh.astype(np.float64) # raw logit scores
y = res.y
scores_neg = yh[y == 0]
scores_pos = yh[y == 1]
print(f"Loaded {SPLIT} split: {len(scores_neg)} negatives, {len(scores_pos)} positives")
print(f"Score range: [{yh.min():.3f}, {yh.max():.3f}]")
# ── Helpers ──────────────────────────────────────────────────────────────────
def ecdf(scores):
"""Return (x, p) for a step-function ECDF."""
s = np.sort(scores)
p = np.arange(1, len(s) + 1) / len(s)
return s, p
def _kde_pdf(scores, x_grid, h):
pdf = FFTKDE(kernel="gaussian", bw=h).fit(scores).evaluate(x_grid)
dx = x_grid[1] - x_grid[0]
cdf = np.cumsum(pdf) * dx
cdf /= cdf[-1]
return cdf
def kde_isj_cdf(scores, x_grid):
"""ISJ bandwidth; Silverman fallback on failure."""
data = scores[:, np.newaxis].astype(np.float64)
try:
h = float(improved_sheather_jones(data))
except Exception:
h = float(silvermans_rule(data))
return _kde_pdf(scores, x_grid, h)
def kde_silverman_cdf(scores, x_grid):
h = float(silvermans_rule(scores[:, np.newaxis].astype(np.float64)))
return _kde_pdf(scores, x_grid, h)
def gauss_cdf(scores, x_grid):
loc, scale = norm.fit(scores)
return norm.cdf(x_grid, loc, scale)
def beta_cdf(scores, x_grid):
"""Fit Beta in sigmoid space; evaluate CDF back on the logit x_grid."""
s01 = expit(scores).clip(1e-6, 1 - 1e-6)
a, b, _, _ = sp_beta.fit(s01, floc=0, fscale=1)
return sp_beta.cdf(expit(x_grid), a, b)
# ── Shared x-grid (covers both classes with a small margin) ──────────────────
margin = 0.5
x_lo = yh.min() - margin
x_hi = yh.max() + margin
x_grid = np.linspace(x_lo, x_hi, 2048)
# ── Compute all curves ───────────────────────────────────────────────────────
neg_ecdf_x, neg_ecdf_p = ecdf(scores_neg)
neg_cdf_kde = kde_isj_cdf(scores_neg, x_grid)
neg_cdf_silverman = kde_silverman_cdf(scores_neg, x_grid)
neg_cdf_gauss = gauss_cdf(scores_neg, x_grid)
neg_cdf_beta = beta_cdf(scores_neg, x_grid)
pos_ecdf_x, pos_ecdf_p = ecdf(scores_pos)
pos_cdf_kde = kde_isj_cdf(scores_pos, x_grid)
pos_cdf_silverman = kde_silverman_cdf(scores_pos, x_grid)
pos_cdf_gauss = gauss_cdf(scores_pos, x_grid)
pos_cdf_beta = beta_cdf(scores_pos, x_grid)
# ── Plot ─────────────────────────────────────────────────────────────────────
fig, (ax0, ax1) = plt.subplots(
2, 1, figsize=(5.5, 4.8), sharex=True,
gridspec_kw={"height_ratios": [1, 1]},
)
fig.subplots_adjust(hspace=0.08)
def _plot_class(ax, color, ecdf_x, ecdf_p,
cdf_isj, cdf_silverman, cdf_gauss, cdf_beta, label_prefix):
ax.step(ecdf_x, ecdf_p,
color=color, lw=1.0, alpha=0.55, where="post",
label="ECDF")
ax.plot(x_grid, cdf_isj,
color=color, lw=1.8, linestyle="-",
label="KDE (ISJ)")
ax.plot(x_grid, cdf_silverman,
color=color, lw=1.4, linestyle="-.",
label="KDE (Silverman)")
ax.plot(x_grid, cdf_gauss,
color=color, lw=1.4, linestyle="--",
label="Gaussian")
ax.plot(x_grid, cdf_beta,
color=color, lw=1.6, linestyle=":",
label="Beta")
ax.set_ylabel("CDF", fontsize=9)
ax.set_yticks([0.0, 0.25, 0.50, 0.75, 1.0])
ax.set_yticklabels(["0", ".25", ".5", ".75", "1"], fontsize=7.5)
ax.tick_params(labelsize=8)
ax.spines[["top", "right"]].set_visible(False)
ax.text(0.97, 0.05, label_prefix,
transform=ax.transAxes, fontsize=8,
ha="right", va="bottom", color=color)
_plot_class(ax0, C_NEG,
neg_ecdf_x, neg_ecdf_p,
neg_cdf_kde, neg_cdf_silverman, neg_cdf_gauss, neg_cdf_beta,
r"Negative class ($y=0$)")
_plot_class(ax1, C_POS,
pos_ecdf_x, pos_ecdf_p,
pos_cdf_kde, pos_cdf_silverman, pos_cdf_gauss, pos_cdf_beta,
r"Positive class ($y=1$)")
ax0.legend(fontsize=7.5, loc="upper left", framealpha=0.88)
ax1.set_xlabel(r"Score $f_\theta(\mathbf{x})$", fontsize=9)
ax0.set_xlim(x_lo, x_hi)
for path_suffix in ("pdf", "png"):
kw = {"bbox_inches": "tight"}
if path_suffix == "png":
kw["dpi"] = 180
fig.savefig(OUT_DIR / f"density_fit_comparison.{path_suffix}", **kw)
print(f"Saved to {OUT_DIR}/density_fit_comparison.{{pdf,png}}")
Binary file not shown.

After

Width:  |  Height:  |  Size: 86 KiB

@@ -0,0 +1,123 @@
"""
Illustrates the gradient locality property of Gaussian-ROLL.
Upper panel: score distributions for both classes with fitted Gaussian PDFs
and operating threshold.
Lower panel: gradient magnitude per sample using the Gaussian gradient formula
(eq:gauss-grad-combined).
Run from the thesis root:
nix develop --command python3 content/method/figures/threshold_weighting_gaussian.py
"""
import numpy as np
from scipy import stats
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
plt.rcParams.update({
"text.usetex": False,
"mathtext.fontset": "cm",
"pdf.fonttype": 42
})
OUT_DIR = "content/method/figures"
np.random.seed(7)
# --- Data ------------------------------------------------------------------
n = 120
alpha = 0.25 # target FPR
neg_scores = np.sort(np.random.normal(-0.6, 0.75, n))
pos_scores = np.sort(np.random.normal(1.1, 0.80, n))
n0 = len(neg_scores)
n1 = len(pos_scores)
# --- Gaussian MLE parameters -----------------------------------------------
mu0, sigma0 = neg_scores.mean(), neg_scores.std()
mu1, sigma1 = pos_scores.mean(), pos_scores.std()
# Threshold: (1-alpha)-quantile of negative class so that FPR = alpha
# (FPR = P(score > tau | y=0) = 1 - F_0(tau), so F_0(tau) = 1-alpha)
tau = stats.norm.ppf(1 - alpha, loc=mu0, scale=sigma0)
# Shared scalar: pdf of positive class evaluated at the threshold
common = (1.0 / (sigma1 * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((tau - mu1) / sigma1) ** 2)
# --- Gaussian ROLL gradients (eq:gauss-grad-combined) ---------------------
# y=1: -(1/n1) * common * (1 + (tau-mu1)*(s-mu1)/sigma1^2)
grad_pos_signed = -(1.0 / n1) * common * (
1 + (tau - mu1) * (pos_scores - mu1) / sigma1 ** 2
)
grad_pos_abs = np.abs(grad_pos_signed)
# y=0: +(1/n0) * common * (1 + (tau-mu0)*(s-mu0)/sigma0^2)
grad_neg_signed = (1.0 / n0) * common * (
1 + (tau - mu0) * (neg_scores - mu0) / sigma0 ** 2
)
grad_neg_abs = np.abs(grad_neg_signed)
# Gradient balance check (should be ~0)
print(f"Gradient balance: pos sum = {grad_pos_signed.sum():.6f}, "
f"neg sum = {grad_neg_signed.sum():.6f}, "
f"total = {(grad_pos_signed.sum() + grad_neg_signed.sum()):.2e}")
# --- Smooth reference curves (continuous gradient functions) ---------------
x_ref = np.linspace(neg_scores.min() - 0.3, pos_scores.max() + 0.3, 400)
smooth_pos_abs = np.abs(-(1.0 / n1) * common * (
1 + (tau - mu1) * (x_ref - mu1) / sigma1 ** 2
))
smooth_neg_abs = np.abs((1.0 / n0) * common * (
1 + (tau - mu0) * (x_ref - mu0) / sigma0 ** 2
))
# --- Plot -----------------------------------------------------------------
C_NEG = "#4477AA"
C_POS = "#CC6633"
C_TAU = "#333333"
fig, (ax1, ax2) = plt.subplots(
2, 1, figsize=(5.5, 4.2), sharex=True,
gridspec_kw={"height_ratios": [1.3, 1.0]}
)
fig.subplots_adjust(hspace=0.06)
# — Top: score distributions with fitted Gaussian PDFs —
ax1.hist(neg_scores, bins=20, density=True, alpha=0.35, color=C_NEG,
label=r"Negative class ($y=0$)")
ax1.hist(pos_scores, bins=20, density=True, alpha=0.35, color=C_POS,
label=r"Positive class ($y=1$)")
ax1.plot(x_ref, stats.norm.pdf(x_ref, mu0, sigma0),
color=C_NEG, linewidth=1.8, label=r"Fitted $\mathcal{N}(\mu_0,\sigma_0^2)$")
ax1.plot(x_ref, stats.norm.pdf(x_ref, mu1, sigma1),
color=C_POS, linewidth=1.8, label=r"Fitted $\mathcal{N}(\mu_1,\sigma_1^2)$")
ax1.axvline(tau, color=C_TAU, linestyle="--", linewidth=1.4,
label=rf"Threshold $\tau$ (FPR $=\alpha={alpha}$)")
ax1.set_ylabel("Density", fontsize=9)
ax1.set_yticks([])
ax1.legend(fontsize=7, loc="upper right", framealpha=0.85)
# — Bottom: gradient magnitudes —
ax2.scatter(neg_scores, grad_neg_abs, color=C_NEG, s=12, zorder=3,
alpha=0.75, label=r"$|\partial\mathcal{L}/\partial f_\theta(\mathbf{x})|$, $y=0$")
ax2.scatter(pos_scores, grad_pos_abs, color=C_POS, s=12, zorder=3,
alpha=0.75, label=r"$|\partial\mathcal{L}/\partial f_\theta(\mathbf{x})|$, $y=1$")
ax2.plot(x_ref, smooth_neg_abs, color=C_NEG, linewidth=1.2, alpha=0.55)
ax2.plot(x_ref, smooth_pos_abs, color=C_POS, linewidth=1.2, alpha=0.55)
ax2.axvline(tau, color=C_TAU, linestyle="--", linewidth=1.4)
ax2.set_ylabel("Gradient magnitude", fontsize=9)
ax2.set_xlabel(r"Score $f_\theta(\mathbf{x})$", fontsize=9)
ax2.set_yticks([])
ax2.legend(fontsize=7.5, loc="upper right", framealpha=0.85)
for ax in (ax1, ax2):
ax.tick_params(labelsize=8)
ax.spines[["top", "right"]].set_visible(False)
out_base = f"{OUT_DIR}/threshold_weighting_gaussian"
fig.savefig(f"{out_base}.pdf", bbox_inches="tight")
fig.savefig(f"{out_base}.png", bbox_inches="tight", dpi=180)
print(f"Saved to {out_base}.{{pdf,png}}")