/******************************************************* * Calculate Incomplete Beta Function Ix(a,b) * * ---------------------------------------------------- * * SAMPLE RUN: * * * * Incomplete Beta Function Ix(a,b) * * * * A = 0.500000 * * B = 5.000000 * * x = 0.200000 * * * * Y = 0.855072 * * * * ---------------------------------------------------- * * Reference: * * "Numerical Recipes, By W.H. Press, B.P. Flannery, * * S.A. Teukolsky and T. Vetterling, Cambridge * * University Press, 1986" [BIBLI 08]. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *******************************************************/ #include #include double a,b,x,y; double GAMMLN(double XX); double BETACF(double A,double B,double X); double BETAI(double A,double B,double X) { // returns the incomplete Beta function Ix(a,b) double BT; if (X < 0.0 || X > 1.0) { printf(" BETAI: Bad argument X (must be 0<=X<=1).\n"); return 0.0; } if (X == 0.0 || X == 1.0) BT=0.0; else BT=exp(GAMMLN(A+B)-GAMMLN(A)-GAMMLN(B)+A*log(X)+B*log(1.0-X)); if (X < (A+1.0)/(A+B+2.0)) return (BT*BETACF(A,B,X)/A); else return (1.0-BT*BETACF(B,A,1.0-X)/B); } double BETACF(double A,double B,double X) { // continued fraction for incomplete Beta function, used by BETAI const int ITMAX = 100; const double EPS = 3e-7; double AM,BM,AZ,QAB,QAP,QAM,BZ,EM,TEM,D,AP,BP,APP,BPP,AOLD; int M; AM=1.0; BM=1.0; AZ=1.0; QAB=A+B; QAP=A+1.0; QAM=A-1.0; BZ=1.0-QAB*X/QAP; for (M=1; M<=ITMAX; M++) { EM=M; TEM=EM+EM; D=EM*(B-M)*X/((QAM+TEM)*(A+TEM)); AP=AZ+D*AM; BP=BZ+D*BM; D=-(A+EM)*(QAB+EM)*X/((A+TEM)*(QAP+TEM)); APP=AP+D*AZ; BPP=BP+D*BZ; AOLD=AZ; AM=AP/BPP; BM=BP/BPP; AZ=APP/BPP; BZ=1.0; if (fabs(AZ-AOLD) < EPS*fabs(AZ)) goto e1; } printf(" BETACF: A or B too big, or ITMAX too small.\n"); return 0.0; e1:return AZ; } double GAMMLN(double XX) { // returns the value Ln(Gamma(XX) for XX>0 double STP,HALF,ONE,FPF,X,TMP,SER; double COF[7]; int J; COF[1]= 76.18009173; COF[2]=-86.50532033; COF[3]=24.01409822; COF[4]=-1.231739516; COF[5]=0.120858003E-2; COF[6]=-0.536382E-5; STP=2.50662827465; HALF=0.5; ONE=1.0; FPF=5.5; X=XX-ONE; TMP=X+FPF; TMP=(X+HALF)*log(TMP)-TMP; SER=ONE; for (J=1; J<=6; J++) { X=X+ONE; SER += COF[J]/X; } return (TMP+log(STP*SER)); } void main() { a=0.50; b=5.0; x=2e-01; y = BETAI(a,b,x); printf("\n Incomplete Beta Function Ix(A,B)\n\n"); printf(" A = %f\n", a); printf(" B = %f\n", b); printf(" x = %f\n\n", x); printf(" Y = %f\n\n", y); } // end of file ibeta.cpp