root / trunk / src / linad99 / vbetai.cpp @ 1131
History  View  Annotate  Download (1002 Bytes)
1 
/**


2 
* $Id: vbetai.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 variable 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 
dvariable betai(const dvariable a, const dvariable b, const dvariable x, int maxit) 
22 
{ 
23 
dvariable 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 
} 