BUG: fix negative samples generated by Wald distribution (#29609)

mmaaz-git · charris · commit b427e83710ab · 2025-10-10T16:22:37.000-06:00
In numpy.random, when generating Wald samples, if the mean and scale
have a large discrepancy, the current implementation suffers from
catastrophic cancellation. An explicit example of this has been added
to test_generator_mt19937.py. The key line implicated in distributions.c:

`X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));`

which has numerical issues when Y &gt;&gt; scale.

I have replaced this with the equivalent rationalized form:
```
d = Y + sqrt(Y) * sqrt(Y + 4 * scale);
X = mean - (2 * mean * Y) / d;
```

And now the test passes.
diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
@@ -845,12 +845,12 @@ double random_noncentral_f(bitgen_t *bitgen_state, double dfnum, double dfden,
 
 double random_wald(bitgen_t *bitgen_state, double mean, double scale) {
   double U, X, Y;
-  double mu_2l;
+  double d;
 
-  mu_2l = mean / (2 * scale);
   Y = random_standard_normal(bitgen_state);
   Y = mean * Y * Y;
-  X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));
+  d = 1 + sqrt(1 + 4 * scale / Y);
+  X = mean * (1 - 2 / d);
   U = next_double(bitgen_state);
   if (U <= mean / (mean + X)) {
     return X;
diff --git a/numpy/random/tests/test_generator_mt19937.py b/numpy/random/tests/test_generator_mt19937.py
@@ -1858,6 +1858,11 @@ def test_wald(self):
                             [2.07093587449261, 0.73073890064369]])
         assert_array_almost_equal(actual, desired, decimal=14)
 
+    def test_wald_nonnegative(self):
+        random = Generator(MT19937(self.seed))
+        samples = random.wald(mean=1e9, scale=2.25, size=1000)
+        assert_(np.all(samples >= 0.0))
+
     def test_weibull(self):
         random = Generator(MT19937(self.seed))
         actual = random.weibull(a=1.23, size=(3, 2))
