Issue #17141: random.vonmisesvariate() no more hangs for large kappas.

This commit is contained in:
Serhiy Storchaka 2013-02-10 19:29:54 +02:00
commit 06e5e730dd
3 changed files with 40 additions and 10 deletions

View File

@ -432,22 +432,20 @@ class Random(_random.Random):
if kappa <= 1e-6: if kappa <= 1e-6:
return TWOPI * random() return TWOPI * random()
a = 1.0 + _sqrt(1.0 + 4.0 * kappa * kappa) s = 0.5 / kappa
b = (a - _sqrt(2.0 * a))/(2.0 * kappa) r = s + _sqrt(1.0 + s * s)
r = (1.0 + b * b)/(2.0 * b)
while 1: while 1:
u1 = random() u1 = random()
z = _cos(_pi * u1) z = _cos(_pi * u1)
f = (1.0 + r * z)/(r + z)
c = kappa * (r - f)
d = z / (r + z)
u2 = random() u2 = random()
if u2 < 1.0 - d * d or u2 <= (1.0 - d) * _exp(d):
if u2 < c * (2.0 - c) or u2 <= c * _exp(1.0 - c):
break break
q = 1.0 / r
f = (q + z) / (1.0 + q * z)
u3 = random() u3 = random()
if u3 > 0.5: if u3 > 0.5:
theta = (mu + _acos(f)) % TWOPI theta = (mu + _acos(f)) % TWOPI

View File

@ -473,6 +473,7 @@ class TestDistributions(unittest.TestCase):
g.random = x[:].pop; g.paretovariate(1.0) g.random = x[:].pop; g.paretovariate(1.0)
g.random = x[:].pop; g.expovariate(1.0) g.random = x[:].pop; g.expovariate(1.0)
g.random = x[:].pop; g.weibullvariate(1.0, 1.0) g.random = x[:].pop; g.weibullvariate(1.0, 1.0)
g.random = x[:].pop; g.vonmisesvariate(1.0, 1.0)
g.random = x[:].pop; g.normalvariate(0.0, 1.0) g.random = x[:].pop; g.normalvariate(0.0, 1.0)
g.random = x[:].pop; g.gauss(0.0, 1.0) g.random = x[:].pop; g.gauss(0.0, 1.0)
g.random = x[:].pop; g.lognormvariate(0.0, 1.0) g.random = x[:].pop; g.lognormvariate(0.0, 1.0)
@ -493,6 +494,7 @@ class TestDistributions(unittest.TestCase):
(g.uniform, (1.0,10.0), (10.0+1.0)/2, (10.0-1.0)**2/12), (g.uniform, (1.0,10.0), (10.0+1.0)/2, (10.0-1.0)**2/12),
(g.triangular, (0.0, 1.0, 1.0/3.0), 4.0/9.0, 7.0/9.0/18.0), (g.triangular, (0.0, 1.0, 1.0/3.0), 4.0/9.0, 7.0/9.0/18.0),
(g.expovariate, (1.5,), 1/1.5, 1/1.5**2), (g.expovariate, (1.5,), 1/1.5, 1/1.5**2),
(g.vonmisesvariate, (1.23, 0), pi, pi**2/3),
(g.paretovariate, (5.0,), 5.0/(5.0-1), (g.paretovariate, (5.0,), 5.0/(5.0-1),
5.0/((5.0-1)**2*(5.0-2))), 5.0/((5.0-1)**2*(5.0-2))),
(g.weibullvariate, (1.0, 3.0), gamma(1+1/3.0), (g.weibullvariate, (1.0, 3.0), gamma(1+1/3.0),
@ -509,8 +511,30 @@ class TestDistributions(unittest.TestCase):
s1 += e s1 += e
s2 += (e - mu) ** 2 s2 += (e - mu) ** 2
N = len(y) N = len(y)
self.assertAlmostEqual(s1/N, mu, places=2) self.assertAlmostEqual(s1/N, mu, places=2,
self.assertAlmostEqual(s2/(N-1), sigmasqrd, places=2) msg='%s%r' % (variate.__name__, args))
self.assertAlmostEqual(s2/(N-1), sigmasqrd, places=2,
msg='%s%r' % (variate.__name__, args))
def test_constant(self):
g = random.Random()
N = 100
for variate, args, expected in [
(g.uniform, (10.0, 10.0), 10.0),
(g.triangular, (10.0, 10.0), 10.0),
#(g.triangular, (10.0, 10.0, 10.0), 10.0),
(g.expovariate, (float('inf'),), 0.0),
(g.vonmisesvariate, (3.0, float('inf')), 3.0),
(g.gauss, (10.0, 0.0), 10.0),
(g.lognormvariate, (0.0, 0.0), 1.0),
(g.lognormvariate, (-float('inf'), 0.0), 0.0),
(g.normalvariate, (10.0, 0.0), 10.0),
(g.paretovariate, (float('inf'),), 1.0),
(g.weibullvariate, (10.0, float('inf')), 10.0),
(g.weibullvariate, (0.0, 10.0), 0.0),
]:
for i in range(N):
self.assertEqual(variate(*args), expected)
def test_von_mises_range(self): def test_von_mises_range(self):
# Issue 17149: von mises variates were not consistently in the # Issue 17149: von mises variates were not consistently in the
@ -526,6 +550,12 @@ class TestDistributions(unittest.TestCase):
msg=("vonmisesvariate({}, {}) produced a result {} out" msg=("vonmisesvariate({}, {}) produced a result {} out"
" of range [0, 2*pi]").format(mu, kappa, sample)) " of range [0, 2*pi]").format(mu, kappa, sample))
def test_von_mises_large_kappa(self):
# Issue #17141: vonmisesvariate() was hang for large kappas
random.vonmisesvariate(0, 1e15)
random.vonmisesvariate(0, 1e100)
class TestModule(unittest.TestCase): class TestModule(unittest.TestCase):
def testMagicConstants(self): def testMagicConstants(self):
self.assertAlmostEqual(random.NV_MAGICCONST, 1.71552776992141) self.assertAlmostEqual(random.NV_MAGICCONST, 1.71552776992141)

View File

@ -250,6 +250,8 @@ Core and Builtins
Library Library
------- -------
- Issue #17141: random.vonmisesvariate() no more hangs for large kappas.
- Issue #17149: Fix random.vonmisesvariate to always return results in - Issue #17149: Fix random.vonmisesvariate to always return results in
[0, 2*math.pi]. [0, 2*math.pi].