Correct InverseGamma's mean and support bounds

beast-base/src/main/java/beast/base/spec/inference/distribution/InverseGamma.java

@@ -125,6 +125,19 @@ public Double inverseCumulativeProbability(double p) {
         return 1.0 / dist.inverseSurvivalProbability(p);
     }
 
+    /**
+     * ScalarDistribution.getMean() delegates to getApacheDistribution().getMean(), i.e. the
+     * internal Gamma(alpha, rate=beta) helper's mean (alpha/beta), not InverseGamma's own mean
+     * of beta/(alpha-1). For alpha ≤ 1 the defining integral diverges; since X > 0 always,
+     * that divergence is to +Infinity (never negative or indeterminate), so +Infinity -- not NaN
+     * -- is returned there, consistent with beta/(alpha-1)'s own limit as alpha -> 1+.
+     */
+    @Override
+    public double getMean() {
+        refresh(); // this make sure distribution parameters are updated if they are sampled during MCMC
+        return alpha > 1 ? beta / (alpha - 1) : Double.POSITIVE_INFINITY;
+    }
+
     @Override
 	public List<Double> sample() {
         if (sampler == null) {
@@ -136,9 +149,30 @@ public List<Double> sample() {
         return List.of(x);
     }
 
+    /**
+     * Support is (0, +Infinity), same as the internal Gamma helper's [0, +Infinity) -- but
+     * stated explicitly here (rather than relying on ScalarDistribution's default, which reads
+     * it off getApacheDistribution()) so it stays correct independent of what that returns.
+     */
+    @Override
+    public Double getLowerBoundOfParameter() {
+        return 0.0;
+    }
+
+    @Override
+    public Double getUpperBoundOfParameter() {
+        return Double.POSITIVE_INFINITY;
+    }
+
+    /**
+     * The internal Gamma helper only exists to draw samples via x = 1/y; it does not represent
+     * InverseGamma's own density/CDF/mean, all of which are overridden above without consulting
+     * it. Returning null (rather than the Gamma object) means any ScalarDistribution method added
+     * in future that is not also overridden here will fail loudly instead of silently returning a
+     * Gamma-distributed answer dressed up as an InverseGamma one.
+     */
     @Override
 	protected Object getApacheDistribution() {
-        refresh(); // this make sure distribution parameters are updated if they are sampled during MCMC
-        return dist;
+        return null;
     }
 } // class InverseGamma

beast-base/src/test/java/beast/base/spec/inference/distribution/InverseGammaTest.java

@@ -5,8 +5,7 @@
 import org.apache.commons.statistics.distribution.GammaDistribution;
 import org.junit.jupiter.api.Test;
 
-import static org.junit.jupiter.api.Assertions.assertEquals;
-import static org.junit.jupiter.api.Assertions.assertNotEquals;
+import static org.junit.jupiter.api.Assertions.*;
 
 /**
  * Simple test for inverse gamma distribution.
@@ -203,4 +202,117 @@ public void testInverseCumulativeProbabilityMatchesInverseGammaNotGamma() {
         assertEquals(0.0, d.inverseCumulativeProbability(0.0), 1e-10);
         assertEquals(Double.POSITIVE_INFINITY, d.inverseCumulativeProbability(1.0));
     }
+
+    /**
+     * ScalarDistribution.getMean() also delegates to getApacheDistribution().getMean(), i.e. the
+     * internal Gamma(alpha, rate=beta) helper's mean (alpha/beta), not InverseGamma's own mean of
+     * beta/(alpha-1), which diverges to +Infinity (never negative or NaN, since X > 0 always)
+     * as alpha -> 1+ and stays there for alpha ≤ 1.
+     */
+    @Test
+    public void testGetMeanMatchesInverseGammaNotGamma() {
+        for (TestData td : tests) {
+            double alpha = td.getShape();
+            double beta = td.getScale();
+
+            InverseGamma d = new InverseGamma();
+            d.initByName("alpha", new RealScalarParam<>(alpha, PositiveReal.INSTANCE),
+                    "beta", new RealScalarParam<>(beta, PositiveReal.INSTANCE));
+
+            GammaDistribution gamma = GammaDistribution.of(alpha, 1.0 / beta);
+
+            assertEquals(beta / (alpha - 1), d.getMean(), 1e-10);
+
+            // it must no longer be the raw Gamma mean (the pre-fix bug)
+            assertNotEquals(gamma.getMean(), d.getMean(), 1e-6);
+        }
+
+        // alpha = 2 is the boundary where the mean is still finite (beta/(alpha-1) = beta) but
+        // the variance (beta^2 / ((alpha-1)^2 (alpha-2))) is not.
+        InverseGamma alphaTwo = new InverseGamma();
+        alphaTwo.initByName("alpha", new RealScalarParam<>(2.0, PositiveReal.INSTANCE),
+                "beta", new RealScalarParam<>(1.0, PositiveReal.INSTANCE));
+        assertEquals(1.0, alphaTwo.getMean(), 1e-10);
+        assertNotEquals(GammaDistribution.of(2.0, 1.0).getMean(), alphaTwo.getMean(), 1e-6);
+
+        // mean diverges to +Infinity, continuously, as alpha approaches 1 from above ...
+        InverseGamma nearOne = new InverseGamma();
+        nearOne.initByName("alpha", new RealScalarParam<>(1.0 + 1e-6, PositiveReal.INSTANCE),
+                "beta", new RealScalarParam<>(2.0, PositiveReal.INSTANCE));
+        assertEquals(2.0 / 1e-6, nearOne.getMean(), 1.0);
+
+        // ... and is +Infinity, not NaN, at and below alpha = 1
+        for (double alpha : new double[]{1.0, 0.9, 0.5}) {
+            InverseGamma d = new InverseGamma();
+            d.initByName("alpha", new RealScalarParam<>(alpha, PositiveReal.INSTANCE),
+                    "beta", new RealScalarParam<>(2.0, PositiveReal.INSTANCE));
+            assertEquals(Double.POSITIVE_INFINITY, d.getMean());
+        }
+    }
+
+    /**
+     * getApacheDistribution() exposes the internal Gamma(alpha, rate=beta) helper, which is only
+     * used for sampling via x = 1/y and does not represent InverseGamma's own density/CDF/mean --
+     * all of which are overridden above without consulting it. It should therefore return null
+     * (per its own documented contract), so any ScalarDistribution method added later that isn't
+     * also overridden here fails loudly rather than silently returning a Gamma-distributed answer.
+     * getLowerBoundOfParameter()/getUpperBoundOfParameter() must therefore also be overridden
+     * directly (not derived from getApacheDistribution()) to keep returning the correct support,
+     * (0, +Infinity), and inverseCumulativeProbability's p<=0/p>=1 boundary handling (which calls
+     * them) must keep working.
+     */
+    @Test
+    public void testGetApacheDistributionIsNull() {
+        InverseGamma d = new InverseGamma();
+        d.initByName("alpha", new RealScalarParam<>(3.0, PositiveReal.INSTANCE),
+                "beta", new RealScalarParam<>(2.0, PositiveReal.INSTANCE));
+
+        assertNull(d.getApacheDistribution());
+
+        assertEquals(0.0, d.getLowerBoundOfParameter(), 1e-10);
+        assertEquals(Double.POSITIVE_INFINITY, d.getUpperBoundOfParameter());
+
+        // inverseCumulativeProbability's boundary handling still works without getApacheDistribution()
+        assertEquals(0.0, d.inverseCumulativeProbability(0.0), 1e-10);
+        assertEquals(Double.POSITIVE_INFINITY, d.inverseCumulativeProbability(1.0));
+    }
+
+    /**
+     * sample() draws y from the internal Gamma(alpha, rate=beta) helper and returns x = 1/y --
+     * this is the original "direct sampler" code path the whole InverseGamma investigation started
+     * from, but it had no test exercising it at all. Checks the empirical mean and empirical CDF
+     * of a large sample against getMean() and the reference CDF values (both independently
+     * verified correct elsewhere in this class), with tolerances set from the distribution's own
+     * (finite, for alpha=3) variance/binomial sampling error so the test isn't flaky.
+     */
+    @Test
+    public void testSampleDrawsFromInverseGamma() {
+        double alpha = 3;
+        double beta = 2;
+
+        InverseGamma d = new InverseGamma();
+        d.initByName("alpha", new RealScalarParam<>(alpha, PositiveReal.INSTANCE),
+                "beta", new RealScalarParam<>(beta, PositiveReal.INSTANCE));
+
+        final int n = 200_000;
+        double sum = 0;
+        int belowOne = 0;
+        for (int i = 0; i < n; i++) {
+            double x = d.sample().getFirst();
+            sum += x;
+            if (x < 1.0) {
+                belowOne++;
+            }
+        }
+        double sampleMean = sum / n;
+        double empiricalCdfAtOne = belowOne / (double) n;
+
+        // Var(X) = beta^2 / ((alpha-1)^2 (alpha-2)) = 1 for alpha=3, beta=2, so SE(mean) ~ 0.0022;
+        // 0.05 is a >20-sigma margin.
+        assertEquals(d.getMean(), sampleMean, 0.05);
+
+        // reference value from TestData for (alpha=3, beta=2) at x=1; SE(proportion) ~ 0.001,
+        // 0.02 is a similarly generous margin.
+        assertEquals(0.67667641618306, empiricalCdfAtOne, 0.02);
+    }
 }