!*********************************************************** !* Simple Linear Regression Minimizing the Sum of Absolute * !* Deviation (SAD) * !* ------------------------------------------------------- * !* SAMPLE RUN: * !* (Find Line Alpha + Beta X, minimizing SAD for next set * !* of data): * !* * !* N = 10 * !* X(1)=1 Y(1)=1.15 * !* X(2)=2 Y(2)=2.22 * !* X(3)=3 Y(3)=2.94 * !* X(4)=4 Y(4)=3.85 * !* X(5)=5 Y(5)=5.10 * !* X(6)=6 Y(6)=5.99 * !* X(7)=7 Y(7)=7.25 * !* X(8)=8 Y(8)=8.04 * !* X(9)=9 Y(9)=9.003 * !* X(10)=10 Y(10)=9.999 * !* * !* Alpha = 0.1667778 * !* Beta = 0.9832222 * !* * !* SAD = 0.8270000 * !* ITER = 3 * !* Error code: 0 * !* * !* ------------------------------------------------------- * !* Ref.: Journal of Applied Statistics (1978) vol.27 no.3. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !*********************************************************** PROGRAM TSIMLP REAL X(10), Y(10), D(10) INTEGER INEXT(10) N=10 DATA X /1.,2.,3.,4.,5.,6.,7.,8.,9.,10./ DATA Y /1.15,2.22,2.94,3.85,5.10,5.99,7.25,8.04,9.003,9.999/ CALL SIMLP(N,X,Y,SAD,ALPHA,BETA,D,ITER,INEXT,IFAULT) print *,' ' print *,' Alpha = ', ALPHA print *,' Beta = ', BETA print *,' ' print *,' SAD = ', SAD print *,' ITER= ', ITER print *,' Error code: ', IFAULT print *,' ' END SUBROUTINE SIMLP(N, X, Y, SAD, ALPHA, BETA, D, ITER, INEXT, IFAULT) ! ALGORITHM AS 132 APPL. STATIST. (1978) VOL.27, NO.3 ! SIMPL: Fit Y = ALPHA + BETA.X + error REAL X(N), Y(N), D(N) INTEGER INEXT(N) DATA ACU/1.0E-06/, BIG/1.0E19/, HALF/0.5/, ZERO/0.0/, ONE/1.0/,TWO/2.0/ ! Initial settings IFAULT = 0 ITER = 0 AHALF = HALF + ACU AONE = AHALF + AHALF ! Determine initial basis D(1) = ZERO Y1 = Y(1) IBAS1 = 1 A1 = X(1) DO 10 I = 2, N IF (ABS(A1 - X(I)) .LT. ACU) GO TO 10 A2 = X(I) IBAS2 = I Y2 = Y(I) GO TO 20 10 CONTINUE IFAULT = 1 RETURN ! Calculate initial beta value 20 DET = ONE / (A2 - A1) AAAA = (A2 * Y1 - A1 * Y2) * DET BBBB = (Y2 - Y1) * DET ! Calculate initial D-vector DO 30 I = 2, N DDD = Y(I) - (AAAA + BBBB * X(I)) D(I) = SIGN(ONE, DDD) 30 CONTINUE TOT1 = ONE TOT2 = X(IBAS2) D(IBAS2) = - ONE DO 40 I = 2, N TOT1 = TOT1 + D(I) TOT2 = TOT2 + D(I) * X(I) 40 CONTINUE T = (A2 * TOT1 - TOT2) * DET IF (ABS(T) .LT. AONE) GO TO 50 DET = - DET GO TO 70 ! Main iterative loop begins 50 T = (TOT2 - A1 * TOT1) * DET IF (ABS(T) .LT. AONE) GO TO 160 IFLAG = 2 IOUT = IBAS2 X(IOUT) = A1 AAA = A1 BBB = A2 GO TO 80 60 T = (TOT2 - A2 * TOT1) * DET IF (ABS(T) .LT. AONE) GO TO 160 70 IFLAG = 1 BBB = A1 AAA = A2 IOUT = IBAS1 80 RHO = SIGN(ONE, T) T = HALF * ABS(T) DET = DET * RHO ! Perform partial sort of ratios INEXT(IBAS1) = IBAS2 RATIO = BIG SUM = AHALF DO 120 I = 1, N DDD = (X(I) - AAA) * DET IF (DDD * D(I) .LE. ACU) GO TO 120 TEST = (Y(I) - AAAA - BBBB * X(I)) / DDD IF (TEST .GE. RATIO) GO TO 120 J = IBAS1 SUM = SUM + ABS(DDD) 90 ISAVE = ABS(INEXT(J)) IF (TEST .GE. D(ISAVE)) GO TO 110 IF (SUM .LT. T) GO TO 100 SUBT = ABS((X(ISAVE) - AAA) * DET) IF (SUM - SUBT .LT. T) GO TO 100 SUM = SUM - SUBT D(ISAVE) = SIGN(1, INEXT(J)) INEXT(J) = INEXT(ISAVE) GO TO 90 100 J = ISAVE ISAVE = ABS(INEXT(J)) IF (TEST .LT. D(ISAVE)) GO TO 100 110 INEXT(I) = INEXT(J) INEXT(J) = SIGN(I, INT(D(I))) D(I) = TEST IF (SUM .LT. T) GO TO 120 IIN = ABS(INEXT(IBAS1)) RATIO = D(IIN) 120 CONTINUE ! Update basic indicators IIN = ABS(INEXT(IBAS1)) J = IIN 130 ISAVE = ABS(INEXT(J)) IF (ISAVE .EQ. IBAS2) GO TO 140 ZZZ = SIGN(1, INEXT(J)) TOT1 = TOT1 - ZZZ - ZZZ TOT2 = TOT2 - TWO * ZZZ * X(ISAVE) D(ISAVE) = - ZZZ J = ISAVE GO TO 130 140 ZZZ = SIGN(1, INEXT(IBAS1)) TOT1 = TOT1 - RHO - ZZZ TOT2 = TOT2 - RHO * BBB - ZZZ * X(IIN) D(IOUT) = - RHO ITER = ITER + 1 IF (IFLAG .EQ. 1) GO TO 150 X(IBAS2) = A2 IBAS2 = IIN D(IBAS2) = - ONE A2 = X(IIN) Y2 = Y(IIN) DET = ONE / (A1 - A2) AAAA = (A1 * Y2 - A2 * Y1) * DET BBBB = (Y1 - Y2) * DET GO TO 60 150 IBAS1 = IIN A1 = X(IIN) D(IBAS1) = ZERO Y1 = Y(IIN) DET = ONE / (A2 - A1) AAAA = (A2 * Y1 - A1 * Y2) * DET BBBB = (Y2 - Y1) * DET GO TO 50 ! Calculate optimal sum of absolute deviations 160 SAD = ZERO DO 170 I = 1, N D(I) = Y(I) - (AAAA + BBBB * X(I)) SAD = SAD + ABS(D(I)) 170 CONTINUE ALPHA = AAAA BETA = BBBB RETURN END !end of file simlp.f90