Add numerical stability test for censored distributions

Test that pm.Censored computes log-probabilities stably at the bounds:
- Right censoring (upper bound): log(1 - CDF) when CDF ≈ 1
- Left censoring (lower bound): log(CDF) when CDF ≈ 0

Uses pm.Censored with Normal(0, 1) at ±40 standard deviations.

Author: maresb

tests/logprob/test_censoring.py

@@ -261,3 +261,49 @@ def test_rounding(rounding_op):
         logprob.eval({xr_vv: test_value}),
         expected_logp,
     )
+
+
+@pytest.mark.parametrize(
+    "censoring_side,bound_value",
+    [
+        ("right", 40.0),  # Far right tail: CDF ≈ 1, need stable log(1-CDF)
+        ("left", -40.0),  # Far left tail: CDF ≈ 0, need stable log(CDF)
+    ],
+)
+def test_censored_logprob_numerical_stability(censoring_side, bound_value):
+    """Test that censored distributions use numerically stable log-probability computations.
+
+    For right-censoring at the upper bound, log(1 - CDF) is computed. When CDF ≈ 1
+    (far right tail), this requires a stable logccdf implementation.
+
+    For left-censoring at the lower bound, log(CDF) is computed. When CDF ≈ 0
+    (far left tail), this requires a stable logcdf implementation.
+
+    This test uses pm.Censored which is the high-level API for censored distributions.
+    """
+    import pymc as pm
+
+    ref_scipy = st.norm(0, 1)
+
+    with pm.Model() as model:
+        normal_dist = pm.Normal.dist(mu=0.0, sigma=1.0)
+        if censoring_side == "right":
+            pm.Censored("y", normal_dist, lower=None, upper=bound_value)
+            expected_logp = ref_scipy.logsf(bound_value)  # log(1 - CDF)
+        else:  # left
+            pm.Censored("y", normal_dist, lower=bound_value, upper=None)
+            expected_logp = ref_scipy.logcdf(bound_value)  # log(CDF)
+
+    # Compile the logp function
+    logp_fn = model.compile_logp()
+
+    # Evaluate at the bound - this is where the log survival/cdf function is used
+    logp_at_bound = logp_fn({"y": bound_value})
+
+    # This should be finite and correct, not -inf
+    assert np.isfinite(logp_at_bound), (
+        f"logp at {censoring_side} bound should be finite, got {logp_at_bound}"
+    )
+    assert np.isclose(logp_at_bound, expected_logp, rtol=1e-6), (
+        f"logp at {censoring_side} bound: got {logp_at_bound}, expected {expected_logp}"
+    )