Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

BUG: stats: make beta and betaprime pdfs work for a larger parameter …

…range
  • Loading branch information...
commit 2ba46dfe7ca00ad574071b44d12aa8327279bfaf 1 parent 47c3c16
Pauli Virtanen pv authored
4 scipy/stats/distributions.py
View
@@ -2254,9 +2254,9 @@ def _rvs(self, a, b):
u2 = gamma.rvs(b,size=self._size)
return (u1 / u2)
def _pdf(self, x, a, b):
- return 1.0/special.beta(a,b)*x**(a-1.0)/(1+x)**(a+b)
+ return np.exp(self._logpdf(x, a, b))
def _logpdf(self, x, a, b):
- return (a-1.0)*log(x) - (a+b)*log(1+x) - log(special.beta(a,b))
+ return special.xlogy(a-1.0, x) - special.xlog1py(a+b, x) - special.betaln(a,b)
def _cdf_skip(self, x, a, b):
# remove for now: special.hyp2f1 is incorrect for large a
x = where(x==1.0, 1.0-1e-6,x)
28 scipy/stats/tests/test_distributions.py
View
@@ -410,6 +410,7 @@ def test_cdf(self):
assert_almost_equal(stats.skellam.cdf(k, mu1, mu2), skcdfR, decimal=5)
+
class TestBeta(TestCase):
def test_logpdf(self):
''' Regression test for Ticket #1326:
@@ -419,7 +420,24 @@ def test_logpdf(self):
assert_almost_equal(logpdf, -0.69314718056)
logpdf = stats.beta.logpdf(0,0.5,1)
assert_almost_equal(logpdf, np.inf)
-
+
+ def test_logpdf_ticket_1866(self):
+ alpha, beta = 267, 1472
+ x = np.array([0.2, 0.5, 0.6])
+ b = stats.beta(alpha, beta)
+ assert_allclose(b.logpdf(x).sum(), -1201.699061824062)
+ assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
+
+
+class TestBetaPrime(TestCase):
+ def test_logpdf(self):
+ alpha, beta = 267, 1472
+ x = np.array([0.2, 0.5, 0.6])
+ b = stats.betaprime(alpha, beta)
+ assert_(np.isfinite(b.logpdf(x)).all())
+ assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
+
+
class TestGamma(TestCase):
def test_pdf(self):
# a few test cases to compare with R
@@ -707,20 +725,20 @@ def test_beta(self):
## (array(6.333333333333333), array(0.055555555555555552))
v = stats.beta.expect(lambda x: (x-19/3.)*(x-19/3.), args=(10,5),
loc=5, scale=2)
- assert_almost_equal(v, 1./18., decimal=14)
+ assert_almost_equal(v, 1./18., decimal=13)
m = stats.beta.expect(lambda x: x, args=(10,5), loc=5., scale=2.)
- assert_almost_equal(m, 19/3., decimal=14)
+ assert_almost_equal(m, 19/3., decimal=13)
ub = stats.beta.ppf(0.95, 10, 10, loc=5, scale=2)
lb = stats.beta.ppf(0.05, 10, 10, loc=5, scale=2)
prob90 = stats.beta.expect(lambda x: 1., args=(10,10), loc=5.,
scale=2.,lb=lb, ub=ub, conditional=False)
- assert_almost_equal(prob90, 0.9, decimal=14)
+ assert_almost_equal(prob90, 0.9, decimal=13)
prob90c = stats.beta.expect(lambda x: 1, args=(10,10), loc=5,
scale=2, lb=lb, ub=ub, conditional=True)
- assert_almost_equal(prob90c, 1., decimal=14)
+ assert_almost_equal(prob90c, 1., decimal=13)
def test_hypergeom(self):
Please sign in to comment.
Something went wrong with that request. Please try again.