scale fix

src/pgm_alternate.c

@@ -87,7 +87,7 @@ static PGM_INLINE float
 bounding_kernel(parameter_t const* pr)
 {
     if (pr->x > pr->t) {
-        static const double a = 0.22579135264472733;  // log(sqrt(pi / 2))
+        static const double a = PGM_LOGPI_2;  // log(pi / 2)
         return expf(pr->h * a + (pr->h - 1.) * pr->logx -
                     PGM_PI2_8 * pr->x - pr->lgammah);
     }

tests/test_polyagamma.py

@@ -219,3 +219,28 @@ def test_log_extreme_value_behaviour():
     assert np.isclose(polyagamma_pdf(1e-3, h=10, z=3, return_log=True), -12473.46649418656)
     assert np.isclose(polyagamma_cdf(1e-16, return_log=True), -1250000000000017.2)
     assert np.isclose(polyagamma_pdf(1e-16, return_log=True), -1249999999999945.8)
+
+
+def _polyagamma_moments(h, z):
+    # closed-form mean and variance of PG(h, z)
+    if z == 0:
+        return h / 4, h / 24
+    t = np.tanh(z / 2)
+    return h * t / (2 * z), h * (np.sinh(z) - z) * (1 - t * t) / (4 * z ** 3)
+
+
+# At these (h, z) the alternate method often accepts from its right-hand
+# proposal. A wrong constant in the right-tail bounding kernel inflated the
+# sample mean by about 1% and the variance by about 5% here.
+@pytest.mark.parametrize("h, z", ((2.0, 1.0), (2.7, 0.0), (3.9, 0.0), (3.9, 1.0)))
+def test_alternate_sample_moments(h, z):
+    rng = np.random.default_rng(1)
+    n = 500_000
+    x = random_polyagamma(h, z, size=n, method="alternate", random_state=rng)
+    mean, var = _polyagamma_moments(h, z)
+    # distance of the sample mean and variance from the exact values, in
+    # standard errors
+    mean_score = (x.mean() - mean) / np.sqrt(x.var() / n)
+    var_score = (x.var() - var) / np.sqrt(np.var((x - x.mean()) ** 2) / n)
+    assert abs(mean_score) < 4, mean_score
+    assert abs(var_score) < 4, var_score