[NumPy-Tickets] [NumPy] #1459: mtrand.beta does not handle small parameters well

[NumPy Trac]

#1459: mtrand.beta does not handle small parameters well
----------------------+-----------------------------------------------------
 Reporter:  JohnReid  |       Owner:  somebody
     Type:  defect    |      Status:  new     
 Priority:  normal    |   Milestone:          
Component:  Other     |     Version:          
 Keywords:            |  
----------------------+-----------------------------------------------------
 I've been using scipy.stats.beta and I've noticed it returns NaN when the
 parameters are small. I've tried using the same parameters in R code and
 it handles them just fine. R also handles parameters up to 17 orders of
 magnitude smaller. Is there any documentation on which parameter ranges
 are acceptable? Can I expect similar results with other distributions?

 Please see my recent post on scipy-user@scipy.org with Josef's follow-up.

 Here's the code I used:

 {{{
 import scipy.stats as S, numpy as N
 from rpy2.robjects import r

 alpha, beta = 0.0710, 0.4222
 for i in xrange(20):
      x_from_scipy = S.beta.rvs(alpha, beta)
      x_from_R = r.rbeta(1, alpha, beta)
      print 'Alpha=%.2e; Beta=%.2e; scipy.stats.beta.rvs=%.2e;
 R.rbeta=%.2e' % (alpha, beta, x_from_scipy, x_from_R[0])
      alpha /= 10.
      beta /= 10.

 }}}

 and the output from it:
 {{{
 Alpha=7.10e-02; Beta=4.22e-01; scipy.stats.beta.rvs=2.75e-11;
 R.rbeta=6.60e-02
 Alpha=7.10e-03; Beta=4.22e-02; scipy.stats.beta.rvs=3.73e-84;
 R.rbeta=4.50e-124
 Alpha=7.10e-04; Beta=4.22e-03; scipy.stats.beta.rvs=1.00e+00;
 R.rbeta=1.00e+00
 Alpha=7.10e-05; Beta=4.22e-04; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-313
 Alpha=7.10e-06; Beta=4.22e-05; scipy.stats.beta.rvs=nan; R.rbeta=1.00e+00
 Alpha=7.10e-07; Beta=4.22e-06; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-315
 Alpha=7.10e-08; Beta=4.22e-07; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-316
 Alpha=7.10e-09; Beta=4.22e-08; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-317
 Alpha=7.10e-10; Beta=4.22e-09; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-318
 Alpha=7.10e-11; Beta=4.22e-10; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-319
 Alpha=7.10e-12; Beta=4.22e-11; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-320
 Alpha=7.10e-13; Beta=4.22e-12; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-321
 Alpha=7.10e-14; Beta=4.22e-13; scipy.stats.beta.rvs=nan; R.rbeta=3.95e-322
 Alpha=7.10e-15; Beta=4.22e-14; scipy.stats.beta.rvs=nan; R.rbeta=1.00e+00
 Alpha=7.10e-16; Beta=4.22e-15; scipy.stats.beta.rvs=nan; R.rbeta=1.00e+00
 Alpha=7.10e-17; Beta=4.22e-16; scipy.stats.beta.rvs=nan; R.rbeta=0.00e+00
 Alpha=7.10e-18; Beta=4.22e-17; scipy.stats.beta.rvs=nan; R.rbeta=1.00e+00
 Alpha=7.10e-19; Beta=4.22e-18; scipy.stats.beta.rvs=nan; R.rbeta=0.00e+00
 Alpha=7.10e-20; Beta=4.22e-19; scipy.stats.beta.rvs=nan; R.rbeta=0.00e+00
 Alpha=7.10e-21; Beta=4.22e-20; scipy.stats.beta.rvs=nan; R.rbeta=0.00e+00
 }}}

-- 
Ticket URL: <http://projects.scipy.org/numpy/ticket/1459>
NumPy <http://projects.scipy.org/numpy>
My example project