8000 Merge pull request #18498 from bashtage/Androp0v-vonmises-fix · numpy/numpy@7d3b555 · GitHub
[go: up one dir, main page]
Skip to content

Commit 7d3b555

Browse files
charrischarris
authored
Merge pull request #18498 from bashtage/Androp0v-vonmises-fix
BUG: Fixed Von Mises distribution for big values of kappa
2 parents d2b969d + aa52959 commit 7d3b555

File tree

7 files changed

+121
-8
lines changed

7 files changed

+121
-8
lines changed

‎numpy/random/include/legacy-distributions.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ extern double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden);
3131
extern double legacy_normal(aug_bitgen_t *aug_state, double loc, double scale);
3232
extern double legacy_standard_gamma(aug_bitgen_t *aug_state, double shape);
3333
extern double legacy_exponential(aug_bitgen_t *aug_state, double scale);
34+
extern double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa);
3435
extern int64_t legacy_random_binomial(bitgen_t *bitgen_state, double p,
3536
int64_t n, binomial_t *binomial);
3637
extern int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n,

‎numpy/random/mtrand.pyx

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,6 @@ cdef extern from "numpy/random/distributions.h":
5151
void random_standard_uniform_fill(bitgen_t* bitgen_state, np.npy_intp cnt, double *out) nogil
5252
int64_t random_positive_int(b 8000 itgen_t *bitgen_state) nogil
5353
double random_uniform(bitgen_t *bitgen_state, double lower, double range) nogil
54-
double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) nogil
5554
double random_laplace(bitgen_t *bitgen_state, double loc, double scale) nogil
5655
double random_gumbel(bitgen_t *bitgen_state, double loc, double scale) nogil
5756
double random_logistic(bitgen_t *bitgen_state, double loc, double scale) nogil
@@ -100,6 +99,7 @@ cdef extern from "include/legacy-distributions.h":
10099
double legacy_f(aug_bitgen_t *aug_state, double dfnum, double dfden) nogil
101100
double legacy_exponential(aug_bitgen_t *aug_state, double scale) nogil
102101
double legacy_power(aug_bitgen_t *state, double a) nogil
102+
double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa) nogil
103103

104104
np.import_array()
105105

@@ -2281,7 +2281,7 @@ cdef class RandomState:
22812281
>>> plt.show()
22822282
22832283
"""
2284-
return cont(&random_vonmises, &self._bitgen, size, self.lock, 2,
2284+
return cont(&legacy_vonmises, &self._bitgen, size, self.lock, 2,
22852285
mu, 'mu', CONS_NONE,
22862286
kappa, 'kappa', CONS_NON_NEGATIVE,
22872287
0.0, '', CONS_NONE, None)

‎numpy/random/src/distributions/distributions.c

Lines changed: 18 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -843,6 +843,7 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
843843
return NPY_NAN;
844844
}
845845
if (kappa < 1e-8) {
846+
/* Use a uniform for very small values of kappa */
846847
return M_PI * (2 * next_double(bitgen_state) - 1);
847848
} else {
848849
/* with double precision rho is zero until 1.4e-8 */
@@ -853,9 +854,23 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
853854
*/
854855
s = (1. / kappa + kappa);
855856
} else {
856-
double r = 1 + sqrt(1 + 4 * kappa * kappa);
857-
double rho = (r - sqrt(2 * r)) / (2 * kappa);
858-
s = (1 + rho * rho) / (2 * rho);
857+
if (kappa <= 1e6) {
858+
/* Path for 1e-5 <= kappa <= 1e6 */
859+
double r = 1 + sqrt(1 + 4 * kappa * kappa);
860+
double rho = (r - sqrt(2 * r)) / (2 * kappa);
861+
s = (1 + rho * rho) / (2 * rho);
862+
} else {
863+
/* Fallback to wrapped normal distribution for kappa > 1e6 */
864+
result = mu + sqrt(1. / kappa) * random_standard_normal(bitgen_state);
865+
/* Ensure result is within bounds */
866+
if (result < -M_PI) {
867+
result += 2*M_PI;
868+
}
869+
if (result > M_PI) {
870+
result -= 2*M_PI;
871+
}
872+
return result;
873+
}
859874
}
860875

861876
while (1) {

‎numpy/random/src/legacy/legacy-distributions.c

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,13 @@
1+
/*
2+
* This file contains generation code for distribution that have been modified
3+
* since Generator was introduced. These are preserved using identical code
4+
* to what was in NumPy 1.16 so that the stream of values generated by
5+
* RandomState is not changed when there are changes that affect Generator.
6+
*
7+
* These functions should not be changed except if they contain code that
8+
* cannot be compiled. They should not be changed for bug fixes, performance
9+
* improvements that can change the values produced, or enhancements to precision.
10+
*/
111
#include "include/legacy-distributions.h"
212

313

@@ -390,3 +400,61 @@ int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) {
390400
binomial_t *binomial) {
391401
return random_multinomial(bitgen_state, n, mnix, pix, d, binomial);
392402
}
403+
404+
double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
405+
double s;
406+
double U, V, W, Y, Z;
407+
double result, mod;
408+
int neg;
409+
if (npy_isnan(kappa)) {
410+
return NPY_NAN;
411+
}
412+
if (kappa < 1e-8) {
413+
return M_PI * (2 * next_double(bitgen_state) - 1);
414+
} else {
415+
/* with double precision rho is zero until 1.4e-8 */
416+
if (kappa < 1e-5) {
417+
/*
418+
* second order taylor expansion around kappa = 0
419+
* precise until relatively large kappas as second order is 0
420+
*/
421+
s = (1. / kappa + kappa);
422+
} else {
423+
/* Path for 1e-5 <= kappa <= 1e6 */
424+
double r = 1 + sqrt(1 + 4 * kappa * kappa);
425+
double rho = (r - sqrt(2 * r)) / (2 * kappa);
426+
s = (1 + rho * rho) / (2 * rho);
427+
}
428+
429+
while (1) {
430+
U = next_double(bitgen_state);
431+
Z = cos(M_PI * U);
432+
W = (1 + s * Z) / (s + Z);
433+
Y = kappa * (s - W);
434+
V = next_double(bitgen_state);
435+
/*
436+
* V==0.0 is ok here since Y >= 0 always leads
437+
* to accept, while Y < 0 always rejects
438+
*/
439+
if ((Y * (2 - Y) - V >= 0) || (log(Y / V) + 1 - Y >= 0)) {
440+
break;
441+
}
442+
}
443+
444+
U = next_double(bitgen_state);
445+
446+
result = acos(W);
447+
if (U < 0.5) {
448+
result = -result;
449+
}
450+
result += mu;
451+
neg = (result < 0);
452+
mod = fabs(result);
453+
mod = (fmod(mod + M_PI, 2 * M_PI) - M_PI);
454+
if (neg) {
455+
mod *= -1;
456+
}
457+
458+
return mod;
459+
}
460+
}

‎numpy/random/tests/test_generator_mt19937.py

Lines changed: 21 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010
assert_warns, assert_no_warnings, assert_array_equal,
1111
assert_array_almost_equal, suppress_warnings)
1212

13-
from numpy.random import Generator, MT19937, SeedSequence
13+
from numpy.random import Generator, MT19937, SeedSequence, RandomState
1414

1515
random = Generator(MT19937())
1616

@@ -1735,6 +1735,26 @@ def test_vonmises_nan(self):
17351735
r = random.vonmises(mu=0., kappa=np.nan)
17361736
assert_(np.isnan(r))
17371737

1738+
@pytest.mark.parametrize("kappa", [1e4, 1e15])
1739+
def test_vonmises_large_kappa(self, kappa):
1740+
random = Generator(MT19937(self.seed))
1741+
rs = RandomState(random.bit_generator)
1742+
state = random.bit_generator.state
1743+
1744+
random_state_vals = rs.vonmises(0, kappa, size=10)
1745+
random.bit_generator.state = state
1746+
gen_vals = random.vonmises(0, kappa, size=10)
1747+
if kappa < 1e6:
1748+
assert_allclose(random_state_vals, gen_vals)
1749+
else:
1750+
assert np.all(random_state_vals != gen_vals)
1751+
1752+
@pytest.mark.parametrize("mu", [-7., -np.pi, -3.1, np.pi, 3.2])
1753+
@pytest.mark.parametrize("kappa", [1e-9, 1e-6, 1, 1e3, 1e15])
1754+
def test_vonmises_large_kappa_range(self, mu, kappa):
1755+
r = random.vonmises(mu, kappa, 50)
1756+
assert_(np.all(r > -np.pi) and np.all(r <= np.pi))
1757+
17381758
def test_wald(self):
17391759
random = Generator(MT19937(self.seed))
17401760
actual = random.wald(mean=1.23, scale=1.54, size=(3, 2))

‎numpy/random/tests/test_generator_mt19937_regressions.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
from numpy.testing import (assert_, assert_array_equal)
22
import numpy as np
33
import pytest
4-
from numpy.random import Generator, MT19937
4+
from numpy.random import Generator, MT19937, RandomState
55

66
mt19937 = Generator(MT19937())
77

88

99
class TestRegression:
1010

11-
def test_VonMises_range(self):
11+
def test_vonmises_range(self):
1212
# Make sure generated random variables are in [-pi, pi].
1313
# Regression test for ticket #986.
1414
for mu in np.linspace(-7., 7., 5):

‎numpy/random/tests/test_randomstate.py

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1238,6 +1238,15 @@ def test_vonmises_small(self):
12381238
r = random.vonmises(mu=0., kappa=1.1e-8, size=10**6)
12391239
assert_(np.isfinite(r).all())
12401240

1241+
def test_vonmises_large(self):
1242+
# guard against changes in RandomState when Generator is fixed
1243+
random.seed(self.seed)
1244+
actual = random.vonmises(mu=0., kappa=1e7, size=3)
1245+
desired = np.array([4.634253748521111e-04,
1246+
3.558873596114509e-04,
1247+
-2.337119622577433e-04])
1248+
assert_array_almost_equal(actual, desired, decimal=8)
1249+
12411250
def test_vonmises_nan(self):
12421251
random.seed(self.seed)
12431252
r = random.vonmises(mu=0., kappa=np.nan)

0 commit comments

Comments
 (0)
0