SUBROUTINE lpest(n, p, x, y, maxit, a, b, sd, r, rate, it, npo, ifault) ! Code converted using TO_F90 by Alan Miller ! Date: 2003-12-17 Time: 15:07:09 ! FORMAL PARAMETERS ! N INTEGER input : the number of points ! P REAL input : p in the Lp norm ! X REAL ARRAY (N) input : the observed values x(i) ! Y REAL ARRAY (N) input : the observed values y(i) ! MAXIT INTEGER input : maximum allowable number of iterations ! A REAL output : the estimate of alfa ! B REAL output : the estimate of beta ! SD REAL output : the Lp norm ! R REAL ARRAY (N) output : the signed residuals ! RATE REAL output : abs(S(k+1)-S(k))/S(k+1), where S(k) is ! the Lp norm of the fit on the kth iteration ! IT INTEGER output : the number of iterations ! NPO INTEGER output : the number of the points on the line ! IFAULT INTEGER output : ! IFAULT = 0 if the routine converged. ! = 1 if return was due to an increase in the norm. ! = 2 if the maximum iterations ! specified was less than 2. ! = 3 if the weighted sample variance of x is 0. ! = 4 if N given less than 2. ! = 5 if convergence has not be achieved within ! the maximum number of iterations specified. ! ALGORITHM AS 110 APPL. STATIST. (1977), VOL.26, NO.1 ! Lp-NORM FIT OF STRAIGHT LINE BY EXTENSION OF SCHLOSSMACHER, ! particularly for 1 <= p <= 2. IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: p REAL, INTENT(IN) :: x(n) REAL, INTENT(IN) :: y(n) INTEGER, INTENT(IN) :: maxit REAL, INTENT(OUT) :: a REAL, INTENT(OUT) :: b REAL, INTENT(OUT) :: sd REAL, INTENT(OUT) :: r(n) REAL, INTENT(OUT) :: rate INTEGER, INTENT(OUT) :: it INTEGER, INTENT(OUT) :: npo INTEGER, INTENT(OUT) :: ifault INTEGER :: i, isw REAL :: absri, a2, b2, div, dx, dxy, eps2, & res, sd2, sw, ssx, spxy, w, xiw, xmean, wp, xi, yi, ymean REAL, PARAMETER :: eps = 1.0E-6 IF(maxit < 2) GO TO 9 IF(n < 2) GO TO 10 ifault = 0 wp = p - 2.0 eps2 = 2.0 * eps sd = 0.0 r(1:n) = 1.0 DO it = 1, maxit npo = 0 ! CALCULATE A AND B BY LEAST SQUARES ON WEIGHTED DATA, ! USING THE HERRAMAN ALGORITHM. ! OMIT OBSERVATIONS WITH SMALL RESIDUALS. sw = 0.0 xmean = 0.0 ymean = 0.0 ssx = 0.0 spxy = 0.0 DO i = 1, n absri = ABS(r(i)) IF(absri <= eps) GO TO 2 w = absri ** wp sw = sw + w div = w / sw xi = x(i) - xmean yi = y(i) - ymean xiw = xi * w dx = xi * xiw dxy = yi * xiw ssx = ssx + dx - dx * div spxy = spxy + dxy - dxy * div xmean = xmean + xi * div ymean = ymean + yi * div CYCLE 2 npo = npo + 1 END DO IF(ssx < eps) GO TO 11 b = spxy / ssx a = ymean - b * xmean ! FORM RESIDUALS AND TEST CONVERGENCE sd2 = 0.0 isw = 0 DO i = 1, n res = y(i) - a - b * x(i) absri = ABS(res) IF(ABS(absri - ABS(r(i))) > eps2) isw = 1 sd2 = sd2 + absri ** p r(i) = res END DO rate = ABS(sd2 - sd) / sd2 IF(isw == 0) RETURN IF(it == 1) GO TO 5 ! TEST FOR INCREASE IN NORM IF(sd2 > sd) GO TO 7 5 sd = sd2 a2 = a b2 = b END DO ! FAILED TO CONVERGE IN MAXIT ITERATIONS ifault = 5 RETURN ! NORM INCREASED, RESTORE A, B, AND R, THEN RETURN 7 ifault = 1 a = a2 b = b2 DO i = 1, n r(i) = y(i) - a - b * x(i) END DO RETURN ! MAXIT SPECIFIED LESS THAN 2 9 ifault = 2 RETURN ! N LESS THAN 2 10 ifault = 4 RETURN ! VARIANCE OF WEIGHTED X IS ZERO 11 ifault = 3 RETURN END SUBROUTINE lpest