ADMB Documentation  11.6rc.3396
 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, double x, int MAXIT);
00012 df1b2variable betacf(const df1b2variable& a,const df1b2variable& b, const df1b2variable& x, int MAXIT);
00013 
00014 //df1b2variable betai(const df1b2variable& a,const df1b2variable& b,
00015  // double x, int maxit=100);
00016 
00028 df1b2variable betai(const df1b2variable & a,const df1b2variable & b,double x, int maxit)
00029 {
00030   df1b2variable bt;
00031 
00032   if (x < 0.0 || x > 1.0) cerr << "Bad x in routine betai" << endl;
00033   if (x == 0.0 || x == 1.0) bt=double(0.0);
00034   else
00035     bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00036   if (x < (value(a)+1.0)/(value(a)+value(b)+2.0))
00037     return bt*betacf(a,b,x,maxit)/a;
00038   else
00039     return 1.0-bt*betacf(b,a,1.0-x,maxit)/b;
00040 }
00041 
00053 df1b2variable pbeta(double x, const df1b2variable & a,const df1b2variable & b, int maxit){
00054   return betai(a,b,x,maxit);
00055 }
00056 
00069 df1b2variable betacf(const df1b2variable& a,const df1b2variable& b,
00070   double x, int MAXIT)
00071 {
00072   int m,m2;
00073   df1b2variable aa,c,d,del,h,qab,qam,qap;
00074 
00075   qab=a+b;
00076   qap=a+double(1.0);
00077   qam=a-double(1.0);
00078   c=double(1.0);
00079   d=double(1.0)-qab*x/qap;
00080   if (fabs(value(d)) < FPMIN) d=FPMIN;
00081   d=double(1.0)/d;
00082   h=d;
00083   for (m=1;m<=MAXIT;m++) {
00084     m2=2*m;
00085     aa=double(m)*(b-double(m))*x/((qam+double(m2))*(a+double(m2)));
00086     d=double(1.0)+aa*d;
00087     if (fabs(value(d)) < FPMIN) d=FPMIN;
00088     c=double(1.0)+aa/c;
00089     if (fabs(value(c)) < FPMIN) c=FPMIN;
00090     d=double(1.0)/d;
00091     h *= d*c;
00092     aa = -(a+double(m))*(qab+double(m))*x/((a+double(m2))*(qap+double(m2)));
00093     d=double(1.0)+aa*d;
00094     if (fabs(value(d)) < FPMIN) d=FPMIN;
00095     c=double(1.0)+aa/c;
00096     if (fabs(value(c)) < FPMIN) c=FPMIN;
00097     d=double(1.0)/d;
00098     del=d*c;
00099     h *= del;
00100     if (fabs(value(del)-double(1.0)) < EPS) break;
00101   }
00102   if (m > MAXIT) cerr << "a or b too big, or MAXIT too small in betacf"
00103          << endl;
00104   return h;
00105 }
00106 
00107 
00108 
00109 
00110 
00111 
00112 df1b2variable pbeta(const df1b2variable & x, const df1b2variable & a, const df1b2variable & b, int maxit)
00113 {
00114   df1b2variable bt;
00115 
00116   if (value(x) < 0.0 || value(x) > 1.0) cerr << "Bad x in routine betai" << endl;
00117   if (value(x) == 0.0 || value(x) == 1.0) bt=double(0.0);
00118   else
00119     bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
00120   if (value(x) < (value(a)+1.0)/(value(a)+value(b)+2.0))
00121     return bt*betacf(a,b,x,maxit)/a;
00122   else
00123     return 1.0-bt*betacf(b,a,1.0-x,maxit)/b;
00124 }
00125 
00126 df1b2variable betacf(const df1b2variable& a,const df1b2variable& b, const df1b2variable& x, int MAXIT)
00127 {
00128   int m,m2;
00129   df1b2variable aa,c,d,del,h,qab,qam,qap;
00130 
00131   qab=a+b;
00132   qap=a+double(1.0);
00133   qam=a-double(1.0);
00134   c=double(1.0);
00135   d=double(1.0)-qab*x/qap;
00136   if (fabs(value(d)) < FPMIN) d=FPMIN;
00137   d=double(1.0)/d;
00138   h=d;
00139   for (m=1;m<=MAXIT;m++) {
00140     m2=2*m;
00141     aa=double(m)*(b-double(m))*x/((qam+double(m2))*(a+double(m2)));
00142     d=double(1.0)+aa*d;
00143     if (fabs(value(d)) < FPMIN) d=FPMIN;
00144     c=double(1.0)+aa/c;
00145     if (fabs(value(c)) < FPMIN) c=FPMIN;
00146     d=double(1.0)/d;
00147     h *= d*c;
00148     aa = -(a+double(m))*(qab+double(m))*x/((a+double(m2))*(qap+double(m2)));
00149     d=double(1.0)+aa*d;
00150     if (fabs(value(d)) < FPMIN) d=FPMIN;
00151     c=double(1.0)+aa/c;
00152     if (fabs(value(c)) < FPMIN) c=FPMIN;
00153     d=double(1.0)/d;
00154     del=d*c;
00155     h *= del;
00156     if (fabs(value(del)-double(1.0)) < EPS) break;
00157   }
00158   if (m > MAXIT) cerr << "a or b too big, or MAXIT too small in betacf"
00159          << endl;
00160   return h;
00161 }
00162 
00163 
00164 
00165 #undef MAXIT
00166 #undef EPS
00167 #undef FPMIN