Statistics
| Revision:

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

History | View | Annotate | Download (1002 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, 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.0-x));
29
  if (x < (a+1.0)/(a+b+2.0))
30
    return bt*betacf(a,b,x,maxit)/a;
31
  else
32
    return 1.0-bt*betacf(b,a,1.0-x,maxit)/b;
33
}