root / trunk / src / linad99 / vbetai.cpp @ 601
History  View  Annotate  Download (1005 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, 
22 
int maxit)

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

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

32 
else

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