'***********************************************************
'* 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<xx<1, the reflection formula can be used first:
'
'     Gamma(1-z) = pi/Gamma(z)/sin(pi*z) = pi*z/Gamma(1+z)/sin(pi*z)
'**********************************************************************
1000 'yy=gammln(xx)
  dim cof(6)
  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#: fpf=5.5#
  x1=xx-one
  tmp=x1+fpf
  tmp=(x1+half)*log(tmp)-tmp
  ser=one
  for j=1 to 6
    x1=x1+one
    ser=ser+cof(j)/x1
  next j
  yy=tmp+log(stp*ser)
return


'***********************************************
'  Returns the incomplete beta function Ix(a,b)
'***********************************************
2000  'zz=betai(a,b,x)
  if x < 0 or x > 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