real function FFINV (alpha, n1, n2, ierr) c Author: Glen Cowan c Date: 5-Jan-1998 c Inverse of the cumulative F distribution. c FFINV uses the following routines (included below) from StatLib: c ALNGAM algorithm AS245 Appl. Statist. (1989) vol. 38, no. 2 c BETAIN algorithm AS63 Appl. Statist. (1973), vol.22, no. 3 c XINBTA algorithm AS109 Appl. Statist. (1977) vol. 26, no. 1 c Modified 26 April 2001 by Glen Cowan: c variable SAE in XINBTA changed from -37D.0 to -308D.0 to avoid c infinite loop (only a problem on digital unix). implicit NONE c arguments integer ierr ! 0 = ok integer n1 integer n2 real alpha c local variables double precision ALNGAM double precision inverse_beta, log_beta double precision one_minus_alpha double precision p, q, r double precision quantile double precision XINBTA c begin FFINV = -1. if (n1 .le. 0 .or. n2 .le. 0) then ierr = 3 return elseif (alpha .lt. 0. .or. alpha .gt. 1.) then ierr = 4 return else ierr = 0 endif one_minus_alpha = 1.0 - alpha p = DFLOAT(n2)/2.d0 q = DFLOAT(n1)/2.d0 r = DFLOAT(n2)/DFLOAT(n1) log_beta = ALNGAM(p,ierr) + ALNGAM(q,ierr) - ALNGAM(p+q,ierr) if ( ierr .ne. 0 ) return inverse_beta = XINBTA(p, q, log_beta, one_minus_alpha, ierr) if ( ierr .ne. 0 ) return quantile = r*(1.d0/inverse_beta - 1.d0) c output is single precision FFINV = quantile return END c----------------------------------------------------------------------- DOUBLE PRECISION FUNCTION ALNGAM(XVALUE, IFAULT) C C ALGORITHM AS245 APPL. STATIST. (1989) VOL. 38, NO. 2 C C Calculation of the logarithm of the gamma function C INTEGER IFAULT DOUBLE PRECISION ALR2PI, FOUR, HALF, ONE, ONEP5, R1(9), R2(9), + R3(9), R4(5), TWELVE, X, X1, X2, XLGE, XLGST, XVALUE, + Y, ZERO C C Coefficients of rational functions C DATA R1/-2.66685 51149 5D0, -2.44387 53423 7D1, + -2.19698 95892 8D1, 1.11667 54126 2D1, + 3.13060 54762 3D0, 6.07771 38777 1D-1, + 1.19400 90572 1D1, 3.14690 11574 9D1, + 1.52346 87407 0D1/ DATA R2/-7.83359 29944 9D1, -1.42046 29668 8D2, + 1.37519 41641 6D2, 7.86994 92415 4D1, + 4.16438 92222 8D0, 4.70668 76606 0D1, + 3.13399 21589 4D2, 2.63505 07472 1D2, + 4.33400 02251 4D1/ DATA R3/-2.12159 57232 3D5, 2.30661 51061 6D5, + 2.74647 64470 5D4, -4.02621 11997 5D4, + -2.29660 72978 0D3, -1.16328 49500 4D5, + -1.46025 93751 1D5, -2.42357 40962 9D4, + -5.70691 00932 4D2/ DATA R4/ 2.79195 31791 8525D-1, 4.91731 76105 05968D-1, + 6.92910 59929 1889D-2, 3.35034 38150 22304D0, + 6.01245 92597 64103D0/ C C Fixed constants C DATA ALR2PI/9.18938 53320 4673D-1/, FOUR/4.D0/, HALF/0.5D0/, + ONE/1.D0/, ONEP5/1.5D0/, TWELVE/12.D0/, ZERO/0.D0/ C C Machine-dependant constants. C A table of values is given at the top of page 399 of the paper. C These values are for the IEEE double-precision format for which C B = 2, t = 53 and U = 1023 in the notation of the paper. C DATA XLGE/5.10D6/, XLGST/1.D+305/ C X = XVALUE ALNGAM = ZERO C C Test for valid function argument C IFAULT = 2 IF (X .GE. XLGST) RETURN IFAULT = 1 IF (X .LE. ZERO) RETURN IFAULT = 0 C C Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined C IF (X .LT. ONEP5) THEN IF (X .LT. HALF) THEN ALNGAM = -LOG(X) Y = X + ONE C C Test whether X < machine epsilon C IF (Y .EQ. ONE) RETURN ELSE ALNGAM = ZERO Y = X X = (X - HALF) - HALF END IF ALNGAM = ALNGAM + X * ((((R1(5)*Y + R1(4))*Y + R1(3))*Y + + R1(2))*Y + R1(1)) / ((((Y + R1(9))*Y + R1(8))*Y + + R1(7))*Y + R1(6)) RETURN END IF C C Calculation for 1.5 <= X < 4.0 C IF (X .LT. FOUR) THEN Y = (X - ONE) - ONE ALNGAM = Y * ((((R2(5)*X + R2(4))*X + R2(3))*X + R2(2))*X + + R2(1)) / ((((X + R2(9))*X + R2(8))*X + R2(7))*X + + R2(6)) RETURN END IF C C Calculation for 4.0 <= X < 12.0 C IF (X .LT. TWELVE) THEN ALNGAM = ((((R3(5)*X + R3(4))*X + R3(3))*X + R3(2))*X + R3(1)) / + ((((X + R3(9))*X + R3(8))*X + R3(7))*X + R3(6)) RETURN END IF C C Calculation for X >= 12.0 C Y = LOG(X) ALNGAM = X * (Y - ONE) - HALF * Y + ALR2PI IF (X .GT. XLGE) RETURN X1 = ONE / X X2 = X1 * X1 ALNGAM = ALNGAM + X1 * ((R4(3)*X2 + R4(2))*X2 + R4(1)) / + ((X2 + R4(5))*X2 + R4(4)) RETURN END c----------------------------------------------------------------------- double precision function xinbta(p,q,beta,alpha,ifault) implicit double precision (a-h,o-z) c c algorithm as 109 appl. statist. (1977), vol.26, no.1 c (replacing algorithm as 64 appl. statist. (1973), c vol.22, no.3) c c Remark AS R83 and the correction in vol40(1) p.236 have been c incorporated in this version. c c Computes inverse of the incomplete beta function c ratio for given positive values of the arguments c p and q, alpha between zero and one. c log of complete beta function, beta, is assumed to be known. c c Auxiliary function required: BETAIN = algorithm AS63 c logical indx c c Define accuracy and initialise. c SAE below is the most negative decimal exponent which does not c cause an underflow; a value of -308 or thereabouts will often be c OK in double precision. c c data acu/1.0d-14/ c data SAE/-37.D0/ ! changed to -308 (GDC, 24 April 2001) data SAE/-308.D0/ data zero/0.0d0/, one/1.0d0/, two/2.0d0/ data three/3.0d0/, four/4.0d0/, five/5.0d0/, six/6.0d0/ c fpu = 10.d0 ** sae xinbta = alpha c c test for admissibility of parameters c ifault = 1 if (p.le.zero .or. q.le.zero) return ifault = 2 if (alpha.lt.zero .or. alpha.gt.one) return ifault = 0 if (alpha.eq.zero .or. alpha.eq.one) return c c change tail if necessary c if (alpha.le.0.5d0) goto 1 a = one-alpha pp = q qq = p indx = .true. goto 2 1 a = alpha pp = p qq = q indx = .false. c c calculate the initial approximation c 2 r = dsqrt(-dlog(a*a)) y = r-(2.30753d0+0.27061d0*r)/(one+(0.99229d0+0.04481d0*r)*r) if(pp.gt.one .and. qq.gt.one) goto 5 r = qq+qq t = one/(9.0d0*qq) t = r*(one-t+y*dsqrt(t))**3 if(t.le.zero) goto 3 t = (four*pp+r-two)/t if(t.le.one) goto 4 xinbta = one-two/(t+one) goto 6 3 xinbta = one-dexp((dlog((one-a)*qq)+beta)/qq) goto 6 4 xinbta = dexp((dlog(a*pp)+beta)/pp) goto 6 5 r = (y*y-three)/six s = one/(pp+pp-one) t = one/(qq+qq-one) h = two/(s+t) w = y*dsqrt(h+r)/h-(t-s)*(r+five/six-two/(three*h)) xinbta = pp/(pp+qq*dexp(w+w)) c c solve for x by a modified newton-raphson method, c using the function betain c 6 r = one-pp t = one-qq yprev = zero sq = one prev = one if(xinbta.lt.0.0001d0) xinbta = 0.0001d0 if(xinbta.gt.0.9999d0) xinbta = 0.9999d0 IEX = MAX(-5.D0/PP**2 - 1.D0/A**2 - 13.D0, SAE) ACU = 10.D0 ** IEX 7 y = betain(xinbta,pp,qq,beta,ifault) if(ifault.eq.0) goto 8 ifault = 3 return 8 continue xin = xinbta y = (y-a)*exp(beta+r*log(xin)+t*log(one-xin)) if(y*yprev.le.zero) prev = max(sq, fpu) g = one 9 adj = g*y sq = adj*adj if(sq.ge.prev) goto 10 tx = xinbta-adj if(tx.ge.zero .and. tx.le.one) goto 11 10 g = g/three goto 9 11 if(prev.le.acu) goto 12 if(y*y.le.acu) goto 12 if(tx.eq.zero .or. tx.eq.one) goto 10 if(tx.eq.xinbta) goto 12 xinbta = tx yprev = y goto 7 12 if (indx) xinbta = one-xinbta return end c----------------------------------------------------------------------- double precision function betain(x, p, q, beta, ifault) implicit double precision (a-h, o-z) c c algorithm as 63 appl. statist. (1973), vol.22, no.3 c c computes incomplete beta function ratio for arguments c x between zero and one, p and q positive. c log of complete beta function, beta, is assumed to be known c logical indx c c define accuracy and initialise c data zero/0.0d0/, one/1.0d0/, acu/0.1d-14/ betain=x c c test for admissibility of arguments c ifault=1 if(p.le.zero .or. q.le.zero) return ifault=2 if(x.lt.zero .or. x.gt.one) return ifault=0 if(x.eq.zero .or. x.eq. one) return c c change tail if necessary and determine s c psq=p+q cx=one-x if(p.ge.psq*x) goto 1 xx=cx cx=x pp=q qq=p indx=.true. goto 2 1 xx=x pp=p qq=q indx=.false. 2 term=one ai=one betain=one ns=qq+cx*psq c c user soper's reduction formulae. c rx=xx/cx 3 temp=qq-ai if(ns.eq.0) rx=xx 4 term=term*temp*rx/(pp+ai) betain=betain+term temp=abs(term) if(temp.le.acu .and. temp.le.acu*betain) goto 5 ai=ai+one ns=ns-1 if(ns.ge.0) goto 3 temp=psq psq=psq+one goto 4 c c calculate result c 5 betain=betain*exp(pp*log(xx)+(qq-one)*log(cx)-beta)/pp if(indx) betain=one-betain return end