Statistics
| Revision:

root / branches / pthreads-737 / src / linad99 / cbetai.cpp @ 765

History | View | Annotate | Download (986 Bytes)

1
/**
2
 * $Id: cbetai.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 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.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
}
34