root / trunk / src / linad99 / cbetai.cpp @ 601
History  View  Annotate  Download (986 Bytes)
1 
/**


2 
* $Id: cbetai.cpp 789 20101005 01:01:09Z johnoel $

3 
*

4 
* Author: David Fournier

5 
* Copyright (c) 2009, 2010 ADMB Foundation

6 
*/

7 
#include <fvar.hpp> 
8 
#include <math.h> 
9  
10 
/** Incomplete beta function for constant objects.

11 
\param a \f$a\f$

12 
\param b \f$b\f$

13 
\param x \f$x\f$

14 
\param maxit Maximum number of iterations for the continued fraction approximation in betacf.

15 
\return Incomplete beta function \f$I_x(a,b)\f$

16 
`

17 
\n\n The implementation of this algorithm was inspired by

18 
"Numerical Recipes in C", 2nd edition,

19 
Press, Teukolsky, Vetterling, Flannery, chapter 2

20 
*/

21 
double betai(const double a,const double b,const double x,int maxit) 
22 
{ 
23 
double bt;

24  
25 
if (x < 0.0  x > 1.0) cerr << "Bad x in routine betai" << endl; 
26 
if (x == 0.0  x == 1.0) bt=0.0; 
27 
else

28 
bt=exp(gammln(a+b)gammln(a)gammln(b)+a*log(x)+b*log(1.0x)); 
29 
if (x < (a+1.0)/(a+b+2.0)) 
30 
return bt*betacf(a,b,x,maxit)/a;

31 
else

32 
return 1.0bt*betacf(b,a,1.0x,maxit)/b; 
33 
} 
34 