diff options
| author | Charles Harris <charlesr.harris@gmail.com> | 2014-05-29 16:30:59 -0600 |
|---|---|---|
| committer | Charles Harris <charlesr.harris@gmail.com> | 2014-05-29 16:30:59 -0600 |
| commit | d1a2f7d92ff00d4e01e3a99124c76b99f561bfc9 (patch) | |
| tree | a490e6475b6cd86519315d80aaac64f15478c4c3 | |
| parent | 3d14321fcec97501a20d3e9ade92abb7c3055c0e (diff) | |
| parent | 69e26e54edd3c83fae61e880a102d1b979e2d67d (diff) | |
| download | numpy-d1a2f7d92ff00d4e01e3a99124c76b99f561bfc9.tar.gz | |
Merge pull request #4751 from juliantaylor/vonmises-fix
BUG: avoid infinite loop for small kappa in vonmises
| -rw-r--r-- | numpy/random/mtrand/distributions.c | 18 | ||||
| -rw-r--r-- | numpy/random/tests/test_random.py | 6 |
2 files changed, 20 insertions, 4 deletions
diff --git a/numpy/random/mtrand/distributions.c b/numpy/random/mtrand/distributions.c index 98f86c064..001e2f6a1 100644 --- a/numpy/random/mtrand/distributions.c +++ b/numpy/random/mtrand/distributions.c @@ -556,7 +556,7 @@ double rk_standard_t(rk_state *state, double df) */ double rk_vonmises(rk_state *state, double mu, double kappa) { - double r, rho, s; + double s; double U, V, W, Y, Z; double result, mod; int neg; @@ -567,9 +567,19 @@ double rk_vonmises(rk_state *state, double mu, double kappa) } else { - r = 1 + sqrt(1 + 4*kappa*kappa); - rho = (r - sqrt(2*r))/(2*kappa); - s = (1 + rho*rho)/(2*rho); + /* with double precision rho is zero until 1.4e-8 */ + if (kappa < 1e-5) { + /* + * second order taylor expansion around kappa = 0 + * precise until relatively large kappas as second order is 0 + */ + s = (1./kappa + kappa); + } + else { + double r = 1 + sqrt(1 + 4*kappa*kappa); + double rho = (r - sqrt(2*r)) / (2*kappa); + s = (1 + rho*rho)/(2*rho); + } while (1) { diff --git a/numpy/random/tests/test_random.py b/numpy/random/tests/test_random.py index ef4e7b127..b64c9d6cd 100644 --- a/numpy/random/tests/test_random.py +++ b/numpy/random/tests/test_random.py @@ -629,6 +629,12 @@ class TestRandomDist(TestCase): [ 1.19153771588353052, 1.83509849681825354]]) np.testing.assert_array_almost_equal(actual, desired, decimal=15) + def test_vonmises_small(self): + # check infinite loop, gh-4720 + np.random.seed(self.seed) + r = np.random.vonmises(mu=0., kappa=1.1e-8, size=10**6) + np.testing.assert_(np.isfinite(r).all()) + def test_wald(self): np.random.seed(self.seed) actual = np.random.wald(mean = 1.23, scale = 1.54, size = (3, 2)) |