'*********************************************************** '* This program calculates the F Distribution Q(F|nu1,nu2) * '* for F, nu1, nu2 given. * '* ------------------------------------------------------- * '* Ref.: "Numerical Recipes, By W.H. Press, B.P. Flannery, * '* S.A. Teukolsky and T. Vetterling, Cambridge University * '* Press, 1986, pages 166-169" [BIBLI 08]. * '* ------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Calculate F Distribution Function Q(F|nu1,nu2) for * '* F, nu1 and nu2 given. * '* * '* Ratio F of dispersion of 1st sample / dispersion of * '* 2nd sample: 1.25 * '* Degree of freedom nu1 of first sample: 25 * '* Degree of freedom nu2 of second sample: 15 * '* * '* Number of iterations: 6 * '* * '* * '* Q = 0.332373 * * '* * '* ------------------------------------------------------- * '* BASIC version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************** defdbl a-h,o-z defint i-n MAXIT=100 EPS=3e-7 FPMIN=1e-30 ' input data: f,xnu1,xnu2 ' we use the formula: Q(F|nu1,nu2) = BetaI(nu2/2,nu1/2,nu2/(nu2+nu1*F)) cls print print " Calculate F Distribution Function Q(F|nu1,nu2) for" print " F, nu1 and nu2 given." print input " Ratio F of dispersion of 1st sample / dispersion of 2nd sample: ",f input " Degree of freedom nu1 of first sample: ",xnu1 input " Degree of freedom nu2 of second sample: ",xnu2 print x=xnu2/(xnu2+xnu1*f) a=xnu2/2#: b=xnu1/2#: gosub 2000 'zz=betai(nu2/2,nu1/2,x) print : print print using " Q = #.######"; zz END '********************************************************************** ' Returns the value ln(Gamma(xx)) for xx>0. Full accuracy is obtained ' for xx > 1. For 0 1 then print " Bad argument x in function BetaI" zz=0: return end if if x=0 or x=1 then bt=0 else 'bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1-x)) xx=a+b: gosub 1000: g1=yy 'g1=gammln(a+b) xx=a: gosub 1000: g2=yy 'g2=gammln(a) xx=b: gosub 1000: g3=yy 'g3=gammln(b) bt=exp(g1-g2-g3+a*log(x)+b*log(1#-x)) end if if x < (a+1#)/(a+b+2#) then 'use continued fraction directly 'zz = bt*betacf(a,b,x)/a aa=a:bb=b:xx=x: gosub 3000: zz=bt*yy/a else 'zz = 1-bt*betacf(b,a,1.-x)/b 'use continued fraction after 'making the symmetry transformation aa=b: bb=a: xx=1#-x: gosub 3000 zz=1#-bt*yy/b end if return '****************************************************************** ' Continued fraction for incomplete beta function (used by BETAI) '****************************************************************** 3000 'yy=betacf(aa,bb,xx) ' int m,m2; ' double a1,c,d,del,h,qab,qam,qap; qab=aa+bb: qap=aa+1#: qam=aa-1# c=1# d=1#-qab*xx/qap if abs(d) < FPMIN then d=FPMIN d=1#/d: h=d for m=1 to MAXIT m2=2*m a1=m*(bb-m)*xx/((qam+m2)*(aa+m2)) d=1#+a1*d if abs(d) < FPMIN then d=FPMIN c=1#+a1/c if abs(c) < FPMIN then c=FPMIN d=1#/d h=h*d*c a1=-(aa+m)*(qab+m)*xx/((qap+m2)*(aa+m2)) d=1#+a1*d if abs(d) < FPMIN then d=FPMIN c=1#+a1/c if abs(c) < FPMIN then c=FPMIN d=1#/d del=d*c h=h*del if abs(del-1#) < EPS then goto 3010 next m 3010 if m>MAXIT then print " a or b too big, or MAXIT too small." else print " Number of iterations: ";m end if yy=h 'return value return 'end of file fdistri.bas