Merge pull request #29926 from charris/backport-29609

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

Author: charris

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;

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))