[SciPy-dev] Numerical Recipes (was tagging 0.7rc1 this weekend?)

[Bruce Southey]

josef.pktd at gmail.com wrote:
>> Hi,
>> Could you also please replace the use of 'betai' with the appropriate
>> call to the t-distribution to the the probabilities?
>> Numerical Recipes uses betai and also it is more understandable to use
>> the actual t-distribution.
>>
>> Bruce
>>
>>     
>
> betai is used in several places in scipy.stats, until now I never
> touched scipy.special calls unless there was a real bug.
>
> There is some duplication in scipy.special, but I have no idea about
> the relative advantages and disadvantages of the different
> implementations. For some functions, I know that they don't work
> correctly over some parameter range. So, until now I'm reluctant to
> change any of the calls to special.
>
> Usually I don't change anything without having some tests first, and
> neither ttest_ind nor ttest_rand have any tests in the test suite and
> I'm running out of time to write tests before 0.70.
>
> but this replacement looks simple enough:
>
>   
>>>> df=5;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
>>>>         
> 0.6382988716409288
> 0.63829887164092902
>   
>>>> df=500;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
>>>>         
> 0.61729505009465102
> 0.61729505009466568
>   
>>>> df=5000;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
>>>>         
> 0.61709708085302761
> 0.61709708085403392
>   
>>>> df=50000;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
>>>>         
> 0.61707727784367095
> 0.61707727785345856
>   
>>>> df=2;t=0.5;stats.betai(0.5*df,0.5,float(df)/(df+t*t));stats.t.sf(t,df)*2
>>>>         
> 0.66666666666666596
> 0.66666666666666652
>
> I will make the change tonight.
>
> Josef
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
>   
The probability density function can be written as an incomplete beta 
function see for example:
http://en.wikipedia.org/wiki/Student%27s_t-distribution

So it is probably just numerical error between division of two gamma 
functions and using the inverse of a beta function. While not an expert 
on this, it looks a like numerical precision of a double since the 
difference starts around 14 decimal place.

Bruce