Statistics
| Revision:

root / trunk / src / linad99 / vbetai.cpp @ 601

History | View | Annotate | Download (1005 Bytes)

1
/**
2
 * $Id: vbetai.cpp 789 2010-10-05 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.0-x));
30
  if (x < (a+1.0)/(a+b+2.0))
31
    return bt*betacf(a,b,x,maxit)/a;
32
  else
33
    return 1.0-bt*betacf(b,a,1.0-x,maxit)/b;
34
}