MODULE rsquared IMPLICIT NONE INTEGER, SAVE :: ifault PUBLIC :: ifault, qr2 CONTAINS !------------------------------------------------------------ ! The following code for calculating quantiles of R^2 (qr2.for) ! is by the author of AS 261 and uses an improved algorithm. FUNCTION qr2(m, size, rho2, p) RESULT(quantile) ! Computes the quantile of the distribution of the square of the ! sample multiple correlation coefficient for given number of ! random variables M, sample size SIZE, square of the population ! multiple correlation coefficient RHO2, and lower tail area P ! Reference: ! Ding, C.G. (1996) `On the computation of the distribution of ! the square of the sample multiple correlation coefficient', ! Comput. Statist. & Data Analysis, vol. 22, 345-350. ! IFAULT is a fault indicator: ! = 1 if there is no convergence after 10 Newton's iterations ! = 2 if any of the input values is illegal ! = 0 otherwise ! No auxiliary algorithm is required ! ELF90 compatible version by Alan Miller ! N.B. The error indicator (ifault) has been put into the module contents ! as functions can only return one result in ELF90. ! Latest revision - 20 May 1997 INTEGER, INTENT(IN) :: m, size REAL, INTENT(IN) :: rho2, p REAL :: quantile ! Local variables REAL :: n, a, b, ab, ga, q0, yp, y, coeff, q, s, t, cdf, pdf, & gb, gab, v, bndcdf, bndpdf, ynew, diff INTEGER :: na, i, nab, nb, iter REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0, & eps = 1.E-06, delta = 1.E-04, rp = 1.772453850905516028 INTEGER, PARAMETER :: itrmax = 10 quantile = p ! Test for admissibility of arguments ifault = 2 IF (m <= 1 .OR. size <= m .OR. rho2 < zero .OR. & rho2 > one .OR. p < zero .OR. p > one) RETURN ifault = 0 IF (p == zero .OR. p == one) RETURN ! Calculate the constants needed for each Newton's iteration a = (m - 1) / two b = (size - m) / two ab = (size - 1) / two IF (MOD(m + 1, 2) == 0) THEN na = a + half ga = one DO i = 1, na ga = ga * i END DO ELSE na = a + one ga = rp DO i = 1, na ga = ga * (i - half) END DO END IF IF (MOD(size - m, 2) == 0) THEN nb = b - half gb = one DO i = 1, nb gb = gb * i END DO ELSE nb = b gb = rp DO i = 1, nb gb = gb * (i - half) END DO END IF IF (MOD(size - 1, 2) == 0) THEN nab = ab - half gab = one DO i = 1, nab gab = gab * i END DO ELSE nab = ab gab = rp DO i = 1, nab gab = gab * (i - half) END DO END IF q0 = (one - rho2) ** ab coeff = gab / ga / gb ! Use 0.5 as a starting value for Newton's iterations y = half ! Perform Newton's iterations DO iter = 1, itrmax ! Evaluate the first terms of the series for CDF (distribution ! function) and PDF (density) n = one yp = one - y t = coeff * y ** a * yp ** b s = a * t / y / yp q = q0 v = q cdf = v * t pdf = q * s ! Check if a + n > (a + b + n)y 80 IF (a + n > (a + b + n) * y) GO TO 90 ! Evaluate the next terms of two series and then the ! partial sums q = q * (a + b + n - one) * rho2 / n v = v + q s = t * (a + b + n - one) / yp t = t * y * (a + b + n - one) / (a + n) cdf = cdf + v * t pdf = pdf + q * s n = n + one GO TO 80 ! Find the error bounds and check for convergence for both ! series 90 bndcdf = t * y * (a + b + n - one) / (a + n - (a + b + n) * y) bndpdf = t * (a + b + n - one) * (one - v) / yp 100 IF (bndcdf <= eps .AND. bndpdf <= eps) GO TO 110 ! Continue to update the terms and then accumulate q = q * (a + b + n - one) * rho2 / n v = v + q IF (bndcdf <= eps) THEN s = s * y * (a + b + n - one) / (a + n - one) pdf = pdf + q * s n = n + one bndpdf = s * y * (a + b + n - one) * (one - v) / (a + n - one) GO TO 100 ELSE IF (bndpdf <= eps) THEN t = t * y * (a + b + n - one) / (a + n) cdf = cdf + v * t n = n + one bndcdf = t * y * (a+b+n-one) / (a+n - (a+b+n)*y) GO TO 100 ELSE s = t * (a + b + n - one) / yp t = t * y * (a + b + n - one) / (a + n) cdf = cdf + v * t pdf = pdf + q * s n = n + one GO TO 90 END IF ! Obtain a new Y and make changes if it is illegal 110 diff = (cdf - p) / pdf ynew = y - diff IF (ynew <= zero) THEN y = y / two ELSE IF (ynew >= one) THEN y = (y + one) / two ELSE y = ynew END IF ! Check for convergence of Newton's iterations IF (ABS(diff) <= delta * y) THEN quantile = y RETURN END IF END DO ifault = 1 RETURN END FUNCTION qr2 END MODULE rsquared PROGRAM main ! This is a driver program that calls QR2 and produces output USE rsquared IMPLICIT NONE INTEGER :: m, size, icode, k REAL :: rho2, p, y 10 WRITE (*, 11) 11 FORMAT (/' ENTER M (>1), N (>M), RHO2 (BETWEEN 0 AND 1), and', & /' P (BETWEEN 0 AND 1) FOR QR2 ==> ') READ (*,*) m, size, rho2, p y = qr2(m, size, rho2, p) icode = ifault + 1 SELECT CASE (icode) CASE (1) WRITE (*, 21) m, size, rho2, p, y 21 FORMAT (/' qr2(', i3, ',', i3, ',', g13.4, ',', g13.4, ') = ', g13.5) CASE (2) WRITE (*, 31) 31 FORMAT (/' NO convergence after 10 Newton ITERATIONS') CASE (3) WRITE (*, 41) 41 FORMAT (/' THE INPUT VALUE IS ILLEGAL!') END SELECT WRITE (*, 61) 61 FORMAT (/ ' ENTER 1 TO CONTINUE OR 0 TO QUIT ==> ') READ (*,*) k IF (k == 1) GO TO 10 STOP END PROGRAM main