Add logistic midpoint characteristics

src/irae_risk/logistic.py

@@ -339,6 +339,126 @@ class LogisticPolyRegression:
 
         return scipy.special.expit(F)
 
+    def get_x50(self, theta):
+        """Return the first predictor value at which the model equals 0.5.
+
+        Since ``expit(F) = 0.5`` exactly when ``F = 0``, ``x50`` is the
+        smallest real root of the polynomial logit
+
+            F(x) = sum_i beta_i(theta) x**i.
+
+        For a monotonic fitted model the root is unique.  Defining "first" as
+        the smallest real root also makes the result unambiguous for an
+        unconstrained polynomial with several 0.5 crossings.
+
+        Parameters
+        ----------
+        theta
+            Model parameter vector of length ``degree + 1``.
+
+        Returns
+        -------
+        float
+            The smallest real solution of ``model(x, theta) = 0.5``.
+
+        Raises
+        ------
+        ValueError
+            If ``theta`` is invalid, the model never reaches 0.5 at a real
+            predictor value, or the model is identically 0.5 and hence has no
+            unique first crossing.
+        """
+        theta = np.asarray(theta, dtype=float)
+        expected_shape = (self.degree + 1,)
+        if theta.shape != expected_shape:
+            raise ValueError(
+                f"theta must have shape {expected_shape}; got {theta.shape}."
+            )
+        if not np.all(np.isfinite(theta)):
+            raise ValueError("theta must contain only finite values.")
+
+        beta = np.asarray(self.get_beta(theta), dtype=float)
+        nonzero = np.flatnonzero(beta != 0.0)
+
+        if nonzero.size == 0:
+            raise ValueError(
+                "x50 is not uniquely defined because the model equals 0.5 "
+                "for every x."
+            )
+
+        polynomial_degree = int(nonzero[-1])
+        if polynomial_degree == 0:
+            raise ValueError("The model does not reach 0.5 for any real x.")
+
+        coefficients = beta[:polynomial_degree + 1]
+
+        # Avoid the unnecessary loss of precision of a general polynomial
+        # root solver in the common linear-logistic case.
+        if polynomial_degree == 1:
+            return float(-coefficients[0] / coefficients[1])
+
+        roots = np.roots(coefficients[::-1])
+        real_roots = []
+        machine_tolerance = 100 * np.finfo(float).eps
+
+        for root in roots:
+            candidate = float(root.real)
+            root_scale = max(1.0, abs(candidate))
+
+            # Repeated real roots may acquire a small imaginary part in a
+            # numerical polynomial-root calculation.  In that case, also
+            # accept the real component when its scaled polynomial residual
+            # is negligible.
+            residual = abs(
+                np.polynomial.polynomial.polyval(candidate, coefficients)
+            )
+            coefficient_scale = np.polynomial.polynomial.polyval(
+                abs(candidate), np.abs(coefficients)
+            )
+            small_imaginary_part = (
+                abs(root.imag) <= machine_tolerance * root_scale
+            )
+            small_residual = residual <= 1e-10 * max(
+                coefficient_scale, np.finfo(float).tiny
+            )
+
+            if small_imaginary_part or small_residual:
+                real_roots.append(candidate)
+
+        if not real_roots:
+            raise ValueError("The model does not reach 0.5 for any real x.")
+
+        return float(min(real_roots))
+
+    def get_s50(self, theta):
+        """Return the probability slope at the model's first 0.5 crossing.
+
+        If ``p(x) = expit(F(x))``, then
+
+            dp/dx = p(x) * (1 - p(x)) * F'(x).
+
+        At ``x50``, ``p(x50) = 0.5``, so the reported midpoint slope is
+        ``F'(x50) / 4``.  The derivative is with respect to the predictor on
+        the scale supplied to the model.
+
+        Parameters
+        ----------
+        theta
+            Model parameter vector of length ``degree + 1``.
+
+        Returns
+        -------
+        float
+            ``d model(x, theta) / dx`` evaluated at ``x = get_x50(theta)``.
+        """
+        x50 = self.get_x50(theta)
+        beta = np.asarray(self.get_beta(theta), dtype=float)
+        derivative_coefficients = np.arange(1, len(beta)) * beta[1:]
+        logit_slope = np.polynomial.polynomial.polyval(
+            x50, derivative_coefficients
+        )
+        return float(logit_slope / 4.0)
+
     """
         Calculate negative log-likelihood function
 

tests/test_logistic.py

@@ -12,6 +12,72 @@ def regression_data():
     return x, y, theta
 
 
+def test_x50_and_s50_for_linear_logistic_model():
+    model = LogisticPolyRegression(degree=1)
+    theta = np.array([-2.0, 0.5])
+
+    x50 = model.get_x50(theta)
+
+    assert x50 == pytest.approx(4.0)
+    assert model.model(np.array([x50]), theta)[0] == pytest.approx(0.5)
+    assert model.get_s50(theta) == pytest.approx(0.5 / 4.0)
+
+
+def test_x50_selects_first_crossing_of_unconstrained_polynomial():
+    model = LogisticPolyRegression(degree=3)
+    # F(x) = (x + 2) (x - 1) (x - 3)
+    theta = np.array([6.0, -5.0, -2.0, 1.0])
+
+    assert model.get_x50(theta) == pytest.approx(-2.0)
+    assert model.get_s50(theta) == pytest.approx(15.0 / 4.0)
+
+
+def test_x50_and_s50_are_invariant_to_monotonic_theta_symmetries():
+    model = LogisticPolyRegression(degree=3, mono=True)
+    theta = np.array([-0.5, 0.4, 1.2, -0.3])
+    equivalent_theta = np.array([-0.5, -0.4, -1.2, 0.3])
+
+    x50 = model.get_x50(theta)
+
+    assert model.model(np.array([x50]), theta)[0] == pytest.approx(0.5)
+    assert model.get_x50(equivalent_theta) == pytest.approx(x50)
+    assert model.get_s50(equivalent_theta) == pytest.approx(
+        model.get_s50(theta)
+    )
+
+
+def test_x50_supports_flat_monotonic_midpoint():
+    model = LogisticPolyRegression(degree=3, mono=True)
+    # F(x) = (x - 1)**3 / 3, so p'(x50) = 0 at x50 = 1.
+    theta = np.array([-1.0 / 3.0, 0.0, 1.0, -1.0])
+
+    x50 = model.get_x50(theta)
+
+    assert x50 == pytest.approx(1.0, abs=1e-5)
+    assert model.model(np.array([x50]), theta)[0] == pytest.approx(0.5)
+    assert model.get_s50(theta) == pytest.approx(0.0, abs=1e-9)
+
+
+def test_x50_rejects_models_without_a_unique_real_midpoint():
+    quadratic = LogisticPolyRegression(degree=2)
+
+    with pytest.raises(ValueError, match="does not reach 0.5"):
+        quadratic.get_x50(np.array([1.0, 0.0, 1.0]))
+
+    with pytest.raises(ValueError, match="not uniquely defined"):
+        quadratic.get_x50(np.zeros(3))
+
+
+def test_x50_validates_theta():
+    model = LogisticPolyRegression(degree=1)
+
+    with pytest.raises(ValueError, match="theta must have shape"):
+        model.get_x50(np.array([1.0]))
+
+    with pytest.raises(ValueError, match="finite"):
+        model.get_x50(np.array([0.0, np.nan]))
+
+
 def test_goodness_of_fit_uses_unpenalized_likelihood_by_default(regression_data):
     x, y, theta = regression_data
     model = LogisticPolyRegression(degree=1, lam=(0.0, 0.5))