ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2bet.cpp
Go to the documentation of this file.
00001 
00007 #include <df1b2fun.h>
00008 //#define EPS double(3.0e-7)
00009 #define EPS double(1.0e-9)
00010 #define FPMIN double(1.0e-30)
00011 df1b2variable betacf(const df1b2variable& a,const df1b2variable& b,
00012   double x, int MAXIT);
00013 
00014 
00015 //df1b2variable betai(const df1b2variable& a,const df1b2variable& b,
00016  // double x, int maxit=100);
00017 
00029 df1b2variable betai(const df1b2variable & a,const df1b2variable & b,double x,
00030   int maxit)
00031 {
00032   df1b2variable bt;
00033 
00034   if (x < 0.0 || x > 1.0) cerr << "Bad x in routine betai" << endl;
00035   if (x == 0.0 || x == 1.0) bt=double(0.0);
00036   else
00037     bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00038   if (x < (value(a)+1.0)/(value(a)+value(b)+2.0))
00039     return bt*betacf(a,b,x,maxit)/a;
00040   else
00041     return 1.0-bt*betacf(b,a,1.0-x,maxit)/b;
00042 }
00043 
00056 df1b2variable betacf(const df1b2variable& a,const df1b2variable& b,
00057   double x, int MAXIT)
00058 {
00059   int m,m2;
00060   df1b2variable aa,c,d,del,h,qab,qam,qap;
00061 
00062   qab=a+b;
00063   qap=a+double(1.0);
00064   qam=a-double(1.0);
00065   c=double(1.0);
00066   d=double(1.0)-qab*x/qap;
00067   if (fabs(value(d)) < FPMIN) d=FPMIN;
00068   d=double(1.0)/d;
00069   h=d;
00070   for (m=1;m<=MAXIT;m++) {
00071     m2=2*m;
00072     aa=double(m)*(b-double(m))*x/((qam+double(m2))*(a+double(m2)));
00073     d=double(1.0)+aa*d;
00074     if (fabs(value(d)) < FPMIN) d=FPMIN;
00075     c=double(1.0)+aa/c;
00076     if (fabs(value(c)) < FPMIN) c=FPMIN;
00077     d=double(1.0)/d;
00078     h *= d*c;
00079     aa = -(a+double(m))*(qab+double(m))*x/((a+double(m2))*(qap+double(m2)));
00080     d=double(1.0)+aa*d;
00081     if (fabs(value(d)) < FPMIN) d=FPMIN;
00082     c=double(1.0)+aa/c;
00083     if (fabs(value(c)) < FPMIN) c=FPMIN;
00084     d=double(1.0)/d;
00085     del=d*c;
00086     h *= del;
00087     if (fabs(value(del)-double(1.0)) < EPS) break;
00088   }
00089   if (m > MAXIT) cerr << "a or b too big, or MAXIT too small in betacf"
00090          << endl;
00091   return h;
00092 }
00093 #undef MAXIT
00094 #undef EPS
00095 #undef FPMIN