fix(#135): make Gaussian profile likelihood scale-safe

Author: marcohost33-maker

src/liouscope/_types.py

@@ -187,6 +187,7 @@ class FitResult:
     n_eff: float
     residual_ar1_rho: float
     success: bool
+    likelihood_degenerate: bool = False
 
 
 @dataclass(frozen=True, slots=True, kw_only=True)

src/liouscope/diagnostics/relaxation.py

@@ -215,6 +215,7 @@ def _fit_with_model(
         n_eff=n_eff,
         residual_ar1_rho=fit.rho_ar1,
         success=fit.success,
+        likelihood_degenerate=fit.likelihood_degenerate,
     )
     return fit_result, fit.params
 

src/liouscope/fitting/aicc.py

@@ -11,10 +11,13 @@
 
 from __future__ import annotations
 
+import math
 from collections.abc import Mapping
 
 import numpy as np
 
+from ..numerics.norms import scaled_log_sum_squares
+
 
 def aicc(log_likelihood: float, k: int, n_eff: float) -> float:
     """Return ``AICc`` for given log-likelihood, parameter count, and N_eff."""
@@ -29,18 +32,48 @@ def gaussian_log_likelihood(
     *,
     sigma: float | None = None,
 ) -> float:
-    """Gaussian log-likelihood for ``y - y_hat``.
+    """Gaussian log-likelihood for ``y - y_hat`` without absolute RSS floors.
 
-    Uses MLE sigma if ``sigma`` is omitted.
+    When ``sigma`` is omitted, evaluate the profile likelihood at the Gaussian
+    MLE ``sigma_hat**2 = RSS / n`` directly in log-RSS space. Exact zero RSS has
+    no finite interior MLE for the positive scale parameter and therefore
+    returns NaN (model-selection likelihood unavailable) instead of inventing
+    an absolute epsilon variance. With an explicitly supplied finite positive
+    ``sigma``, zero residuals remain a valid finite likelihood.
     """
     residuals = np.asarray(residuals, dtype=float)
     n = residuals.size
-    rss = float(np.dot(residuals, residuals))
+    if n == 0:
+        return float("nan")
+
+    log_rss = scaled_log_sum_squares(residuals)
+    if math.isnan(log_rss):
+        return float("nan")
+
+    log_2pi = math.log(2.0 * math.pi)
     if sigma is None:
-        sigma_sq = max(rss / n, 1.0e-30)
+        if log_rss == float("-inf"):
+            return float("nan")
+        if log_rss == float("inf"):
+            return float("-inf")
+        return float(-0.5 * n * (log_2pi + 1.0 + log_rss - math.log(n)))
+
+    sigma = float(sigma)
+    if not math.isfinite(sigma) or sigma <= 0.0:
+        return float("nan")
+    log_sigma = math.log(sigma)
+    if log_rss == float("-inf"):
+        standardised_rss = 0.0
+    elif log_rss == float("inf"):
+        return float("-inf")
     else:
-        sigma_sq = sigma * sigma
-    return float(-0.5 * n * (np.log(2.0 * np.pi * sigma_sq) + rss / (n * sigma_sq)))
+        log_standardised_rss = log_rss - 2.0 * log_sigma
+        if log_standardised_rss > math.log(np.finfo(float).max):
+            return float("-inf")
+        standardised_rss = math.exp(log_standardised_rss)
+    return float(
+        -0.5 * n * log_2pi - n * log_sigma - 0.5 * standardised_rss
+    )
 
 
 def choose_model(aiccs: Mapping[str, float]) -> str:

src/liouscope/fitting/bootstrap.py

@@ -65,6 +65,8 @@ def parametric_bootstrap(
             + (f" (saturated: {', '.join(base.saturated)})" if base.saturated else "")
             + (" (the curve carries no resolvable variation, issue #123)"
                if base.degenerate else "")
+            + (" (the residual likelihood scale is degenerate, issue #135)"
+               if base.likelihood_degenerate else "")
             + "; a bootstrap around a non-estimate has no meaning"
         )
     theta_hat = base.params

src/liouscope/fitting/gls.py

@@ -12,13 +12,15 @@
 
 from __future__ import annotations
 
+import math
 import warnings
 from collections.abc import Callable
 from dataclasses import dataclass
 
 import numpy as np
 from scipy.optimize import least_squares
 
+from ..numerics.norms import scaled_log_sum_squares
 from .aicc import gaussian_log_likelihood
 from .models import saturation_watch
 from .neff import _AR1_SMALL_N, ar1_correlation_corrected
@@ -41,6 +43,11 @@ class GLSFitOutput:
     #: reports a fit that ran and ended on a magnitude plateau; here there was
     #: nothing to fit. Implies ``success`` is False and ``params`` is NaN.
     degenerate: bool = False
+    #: True when the residual Gaussian scale has no finite usable MLE for model
+    #: selection (issue #135). Distinct from ``degenerate`` above: the curve
+    #: may carry variation and the optimiser may have run, but exact-zero RSS
+    #: (or an unrepresentable positive MLE scale) cannot support AICc/CI claims.
+    likelihood_degenerate: bool = False
 
 
 def _whiten(y: np.ndarray, rho: float) -> np.ndarray:
@@ -192,14 +199,55 @@ def residual(params: np.ndarray, rho_local: float = rho) -> np.ndarray:
         )
     whitened = _whiten(residuals_raw, rho)
     n = whitened.size
-    sigma = float(np.sqrt(max(np.dot(whitened, whitened) / max(n, 1), 1.0e-30)))
+    log_rss = scaled_log_sum_squares(whitened)
+    if log_rss == float("-inf") or not math.isfinite(log_rss):
+        warnings.warn(
+            "fit_gls_ar1: residual Gaussian scale has no finite positive MLE "
+            "for model selection; likelihood/AICc/CI evidence is unavailable "
+            "(issue #135).",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+        return GLSFitOutput(
+            params=p,
+            residuals=residuals_raw,
+            rho_ar1=rho,
+            sigma=float("nan"),
+            log_likelihood=float("nan"),
+            success=False,
+            saturated=tuple(sorted(fired)),
+            likelihood_degenerate=True,
+        )
+
+    log_sigma = 0.5 * (log_rss - math.log(n))
+    try:
+        sigma = float(math.exp(log_sigma))
+    except OverflowError:
+        sigma = float("inf")
+    if not math.isfinite(sigma) or sigma <= 0.0:
+        warnings.warn(
+            "fit_gls_ar1: positive residual MLE scale is not representable as "
+            "float64; likelihood/AICc/CI evidence is unavailable (issue #135).",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+        return GLSFitOutput(
+            params=p,
+            residuals=residuals_raw,
+            rho_ar1=rho,
+            sigma=float("nan"),
+            log_likelihood=float("nan"),
+            success=False,
+            saturated=tuple(sorted(fired)),
+            likelihood_degenerate=True,
+        )
+
     # Prais-Winsten exact AR(1) likelihood: _whiten keeps observation 0
     # (scaled by sqrt(1-rho^2)), so the transform has log-Jacobian
-    # 0.5*log(1-rho^2). Omitting it makes the reported value a hybrid of the
-    # exact and conditional likelihoods and biases cross-model AICc (each
-    # model fits its own rho) toward under-fitting high-rho models.
+    # 0.5*log(1-rho^2). The profile likelihood is evaluated directly from
+    # log(RSS), not by squaring ``sigma`` or materialising RSS.
     jac = 0.5 * float(np.log(max(1.0 - rho * rho, 1.0e-12)))
-    log_lik = gaussian_log_likelihood(whitened, sigma=sigma) + jac
+    log_lik = gaussian_log_likelihood(whitened) + jac
     return GLSFitOutput(
         params=p,
         residuals=residuals_raw,

src/liouscope/numerics/norms.py

@@ -73,6 +73,40 @@ def scaled_euclidean_norm(values: np.ndarray) -> float:
         return float(np.ldexp(scaled_norm, exponent))
 
 
+def scaled_log_sum_squares(values: np.ndarray) -> float:
+    """Return ``log(sum(abs(values)**2))`` without spurious under/overflow.
+
+    The return value is ``-inf`` for an exact all-zero input, ``nan`` when any
+    component is NaN, and ``inf`` when any component is infinite. For finite,
+    non-zero float64 input the logarithm remains finite even when the true sum
+    of squares (or its square root) is outside the representable float64 range.
+
+    This is the likelihood-facing companion to :func:`scaled_euclidean_norm`:
+    both use the same exact power-of-two scaling contract, but this function
+    never reconstructs RSS in ordinary floating-point units.
+    """
+    arr = np.asarray(values)
+    if arr.size == 0:
+        return float("-inf")
+
+    real = np.asarray(np.real(arr), dtype=float)
+    imag = np.asarray(np.imag(arr), dtype=float)
+    if np.any(np.isnan(real)) or np.any(np.isnan(imag)):
+        return float("nan")
+    if np.any(np.isinf(real)) or np.any(np.isinf(imag)):
+        return float("inf")
+
+    scaled = _finite_component_scale(arr)
+    if scaled is None:
+        return float("-inf")
+    scaled_real, scaled_imag, exponent = scaled
+    sumsq = float(
+        np.sum(scaled_real * scaled_real, dtype=float)
+        + np.sum(scaled_imag * scaled_imag, dtype=float)
+    )
+    return float(math.log(sumsq) + 2.0 * exponent * math.log(2.0))
+
+
 def scaled_cancellation_ratio(values: np.ndarray) -> float:
     """Return ``abs(sum(values)) / sum(abs(values))`` scale-safely.
 

tests/test_issue135_likelihood_scale.py

@@ -0,0 +1,113 @@
+"""Issue #135: scale-safe Gaussian likelihood and degenerate RSS semantics."""
+
+from __future__ import annotations
+
+import math
+
+import numpy as np
+import pytest
+
+from liouscope.diagnostics import relaxation as relaxation_mod
+from liouscope.fitting.aicc import aicc, choose_model, gaussian_log_likelihood
+from liouscope.fitting.bootstrap import parametric_bootstrap
+from liouscope.fitting.gls import GLSFitOutput, fit_gls_ar1
+from liouscope.fitting.models import M0
+from liouscope.numerics.norms import scaled_log_sum_squares
+
+
+def test_scaled_log_rss_spans_underflow_and_overflow_regimes():
+    tiny = np.array([1.0e-320, -2.0e-320])
+    huge = np.array([1.0e308, -1.0e308])
+    log_tiny = scaled_log_sum_squares(tiny)
+    log_huge = scaled_log_sum_squares(huge)
+    assert np.isfinite(log_tiny)
+    assert np.isfinite(log_huge)
+    assert log_tiny == pytest.approx(math.log(5.0) + 2.0 * math.log(1.0e-320), rel=2e-5)
+    assert log_huge == pytest.approx(math.log(2.0) + 2.0 * math.log(1.0e308), rel=1e-14)
+    assert scaled_log_sum_squares(np.zeros(4)) == float("-inf")
+
+
+def test_profile_loglikelihood_is_scale_covariant_and_delta_invariant():
+    x = np.linspace(0.2, 1.2, 64)
+    r0 = 0.7 + x + 0.03 * np.sin(5.0 * x)
+    r1 = 0.9 + 1.1 * x - 0.02 * np.cos(3.0 * x)
+    reference_delta = None
+    reference_winner = None
+    for scale in (1.0e-150, 1.0e-40, 1.0, 1.0e40, 1.0e150):
+        ll0 = gaussian_log_likelihood(scale * r0)
+        ll1 = gaussian_log_likelihood(scale * r1)
+        assert np.isfinite(ll0) and np.isfinite(ll1)
+        delta = ll1 - ll0
+        scores = {
+            "M0": aicc(ll0, k=2, n_eff=64.0),
+            "M1": aicc(ll1, k=3, n_eff=64.0),
+        }
+        winner = choose_model(scores)
+        if reference_delta is None:
+            reference_delta = delta
+            reference_winner = winner
+        else:
+            assert delta == pytest.approx(reference_delta, rel=1e-11, abs=1e-10)
+            assert winner == reference_winner
+
+
+def test_profile_loglikelihood_handles_true_rss_above_float_range():
+    ll = gaussian_log_likelihood(np.array([1.0e308, -1.0e308]))
+    assert np.isfinite(ll)
+
+
+def test_zero_rss_unknown_sigma_is_unavailable_but_known_sigma_is_valid():
+    residuals = np.zeros(8)
+    assert np.isnan(gaussian_log_likelihood(residuals))
+    ll = gaussian_log_likelihood(residuals, sigma=2.0)
+    expected = -0.5 * residuals.size * math.log(2.0 * math.pi * 4.0)
+    assert ll == pytest.approx(expected, rel=1e-14, abs=1e-14)
+
+
+def test_gls_exact_fit_marks_likelihood_degenerate_and_bootstrap_refuses():
+    t = np.linspace(0.0, 4.0, 48)
+    p0 = np.array([1.25, 0.6])
+    y = M0(t, p0)
+    with pytest.warns(RuntimeWarning, match="likelihood/AICc/CI evidence is unavailable"):
+        fit = fit_gls_ar1(M0, t, y, p0, n_iters=1)
+    assert not fit.success
+    assert fit.likelihood_degenerate
+    assert not fit.degenerate
+    assert np.isnan(fit.sigma)
+    assert np.isnan(fit.log_likelihood)
+    # The repository treats unexpected warnings as errors. Here the warning is
+    # part of the intended public contract: the base fit first reports why its
+    # likelihood evidence is unusable, then bootstrap refuses the non-estimate.
+    with pytest.warns(RuntimeWarning, match="likelihood/AICc/CI evidence is unavailable"):
+        with pytest.raises(RuntimeError, match="likelihood scale is degenerate"):
+            parametric_bootstrap(M0, t, y, p0, B=4)
+
+
+def test_ordinary_noisy_gls_positive_control_is_not_likelihood_degenerate(rng):
+    t = np.linspace(0.0, 4.0, 80)
+    p = np.array([1.25, 0.6])
+    y = M0(t, p) + 1.0e-3 * rng.standard_normal(t.size)
+    fit = fit_gls_ar1(M0, t, y, p, n_iters=1)
+    assert fit.success
+    assert not fit.likelihood_degenerate
+    assert np.isfinite(fit.sigma) and fit.sigma > 0.0
+    assert np.isfinite(fit.log_likelihood)
+
+
+def test_likelihood_degenerate_state_reaches_fitresult_and_is_nonselectable(monkeypatch):
+    t = np.linspace(0.0, 1.0, 16)
+    y = np.exp(-t)
+    fake = GLSFitOutput(
+        params=np.array([1.0, 1.0]),
+        residuals=np.zeros_like(y),
+        rho_ar1=0.0,
+        sigma=float("nan"),
+        log_likelihood=float("nan"),
+        success=False,
+        likelihood_degenerate=True,
+    )
+    monkeypatch.setattr(relaxation_mod, "fit_gls_ar1", lambda *args, **kwargs: fake)
+    fit_result, _ = relaxation_mod._fit_with_model("M0", t, y)
+    assert not fit_result.success
+    assert fit_result.likelihood_degenerate
+    assert np.isinf(fit_result.aicc)