Skip to content
This repository

BUG: stats: make beta and betaprime pdfs work for a larger parameter range #473

Closed
wants to merge 1 commit into from

2 participants

Pauli Virtanen Per A. Brodtkorb
Pauli Virtanen
Owner
pv commented March 16, 2013

Gets rid of spurious 0/0 occurrence in the formula.

Closes Trac #1866

Per A. Brodtkorb
pbrod commented March 16, 2013

Should also fix the 0*log(0) situations: e.g.
beta.logpdf= special.xlogy(a-1, x) + special.xlog1py(b-1, -x) - special.betaln(a,b)

similiar change for betaprime

See #467

Pauli Virtanen
Owner
pv commented March 21, 2013

Went in dbaaad3

Pauli Virtanen pv closed this March 21, 2013
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Showing 1 unique commit by 1 author.

Mar 16, 2013
Pauli Virtanen BUG: stats: make beta and betaprime pdfs work for a larger parameter …
…range

Fixes Trac #1866
3d091d0
This page is out of date. Refresh to see the latest.
10  scipy/stats/distributions.py
@@ -2186,12 +2186,10 @@ class beta_gen(rv_continuous):
2186 2186
     def _rvs(self, a, b):
2187 2187
         return mtrand.beta(a,b,self._size)
2188 2188
     def _pdf(self, x, a, b):
2189  
-        Px = (1.0-x)**(b-1.0) * x**(a-1.0)
2190  
-        Px /= special.beta(a,b)
2191  
-        return Px
  2189
+        return np.exp(self._logpdf(x, a, b))
2192 2190
     def _logpdf(self, x, a, b):
2193 2191
         lPx = (b-1.0)*log(1.0-x) + (a-1.0)*log(x)
2194  
-        lPx -= log(special.beta(a,b))
  2192
+        lPx -= special.betaln(a,b)
2195 2193
         return lPx
2196 2194
     def _cdf(self, x, a, b):
2197 2195
         return special.btdtr(a,b,x)
@@ -2256,9 +2254,9 @@ def _rvs(self, a, b):
2256 2254
         u2 = gamma.rvs(b,size=self._size)
2257 2255
         return (u1 / u2)
2258 2256
     def _pdf(self, x, a, b):
2259  
-        return 1.0/special.beta(a,b)*x**(a-1.0)/(1+x)**(a+b)
  2257
+        return np.exp(self._logpdf(x, a, b))
2260 2258
     def _logpdf(self, x, a, b):
2261  
-        return (a-1.0)*log(x) - (a+b)*log(1+x) - log(special.beta(a,b))
  2259
+        return (a-1.0)*log(x) - (a+b)*log(1+x) - special.betaln(a,b)
2262 2260
     def _cdf_skip(self, x, a, b):
2263 2261
         # remove for now: special.hyp2f1 is incorrect for large a
2264 2262
         x = where(x==1.0, 1.0-1e-6,x)
22  scipy/stats/tests/test_distributions.py
@@ -693,20 +693,20 @@ def test_beta(self):
693 693
 ##        (array(6.333333333333333), array(0.055555555555555552))
694 694
         v = stats.beta.expect(lambda x: (x-19/3.)*(x-19/3.), args=(10,5),
695 695
                               loc=5, scale=2)
696  
-        assert_almost_equal(v, 1./18., decimal=14)
  696
+        assert_almost_equal(v, 1./18., decimal=13)
697 697
 
698 698
         m = stats.beta.expect(lambda x: x, args=(10,5), loc=5., scale=2.)
699  
-        assert_almost_equal(m, 19/3., decimal=14)
  699
+        assert_almost_equal(m, 19/3., decimal=13)
700 700
 
701 701
         ub = stats.beta.ppf(0.95, 10, 10, loc=5, scale=2)
702 702
         lb = stats.beta.ppf(0.05, 10, 10, loc=5, scale=2)
703 703
         prob90 = stats.beta.expect(lambda x: 1., args=(10,10), loc=5.,
704 704
                                    scale=2.,lb=lb, ub=ub, conditional=False)
705  
-        assert_almost_equal(prob90, 0.9, decimal=14)
  705
+        assert_almost_equal(prob90, 0.9, decimal=13)
706 706
 
707 707
         prob90c = stats.beta.expect(lambda x: 1, args=(10,10), loc=5,
708 708
                                     scale=2, lb=lb, ub=ub, conditional=True)
709  
-        assert_almost_equal(prob90c, 1., decimal=14)
  709
+        assert_almost_equal(prob90c, 1., decimal=13)
710 710
 
711 711
 
712 712
     def test_hypergeom(self):
@@ -1031,5 +1031,19 @@ def test_ncx2_tails_ticket_955():
1031 1031
     b = stats.ncx2.veccdf(np.arange(20, 25, 0.2), 2, 1.07458615e+02)
1032 1032
     assert_allclose(a, b, rtol=1e-3, atol=0)
1033 1033
 
  1034
+def test_beta_logpdf_ticket_1866():
  1035
+    alpha, beta = 267, 1472
  1036
+    x = np.array([0.2, 0.5, 0.6])
  1037
+    b = stats.beta(alpha, beta)
  1038
+    assert_allclose(b.logpdf(x).sum(), -1201.699061824062)
  1039
+    assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
  1040
+
  1041
+def test_betaprime_logpdf():
  1042
+    alpha, beta = 267, 1472
  1043
+    x = np.array([0.2, 0.5, 0.6])
  1044
+    b = stats.betaprime(alpha, beta)
  1045
+    assert_(np.isfinite(b.logpdf(x)).all())
  1046
+    assert_allclose(b.pdf(x), np.exp(b.logpdf(x)))
  1047
+
1034 1048
 if __name__ == "__main__":
1035 1049
     run_module_suite()
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.