Go to the documentation of this file.00001
00007 #include <df1b2fun.h>
00008
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
00016
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
1.8.0
