'*********************************************************** '* 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 = .1667777962154812 * '* Beta = .9832221799426608 * '* * '* SAD = .8269999027252202 * '* ITER = 3 * '* Error code: 0 * '* * '* ------------------------------------------------------- * '* Ref.: Journal of Applied Statistics (1978) vol.27 no.3. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '*********************************************************** 'PROGRAM TSIMLP DEFDBL A-H, O-Z DEFINT I-N ACU = .000001# BIG = 1D+19 HALF = .5 ZERO = 0# ONE = 1# TWO = 2# ' N: Number of points ' X(N), Y(N): Input Data Set ' D(N): Utility vector ' INEXT(N): idem ' SAD: Sum of Absolute Deviation ' Alpha,Beta: Coefficients of linear regression ' Iter: Number of iterations ' Ifault: Error code (must be zero) N = 10 DIM X(N), Y(N), D(N), INEXT(N) FOR I = 1 TO N X(I) = 1# * I NEXT I Y(1) = 1.15: Y(2) = 2.22: Y(3) = 2.94: Y(4) = 3.85: Y(5) = 5.1 Y(6) = 5.99: Y(7) = 7.25: Y(8) = 8.04: Y(9) = 9.003: Y(10) = 9.999 GOSUB 1000 'call SIMLP(N,X,Y,SAD,ALPHA,BETA,D,ITER,INEXT,IFAULT) CLS PRINT PRINT " Alpha = "; ALPHA PRINT " Beta = "; beta PRINT PRINT " SAD = "; SAD PRINT " ITER = "; ITER PRINT " Error code: "; IFAULT PRINT END 500 'Function ISign(ia,ib) IF ib < 0 THEN ISign = -ABS(ia) ELSE ISign = ABS(ia) END IF RETURN 600 'Function Sign(a,b) IF b < 0# THEN Sign = -ABS(a) ELSE Sign = ABS(a) END IF RETURN 1000 '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 ' A1,A2,AAAA,BBBB,AHALF,AONE,DDD,DET,TOT1,TOT2,Y1,Y2: Double ' AAA,BBB,RATIO,RHO,SUBT,SUM,T,TEST,ZZZ: Double ' I,IBAS1,IBAS2,IFLAG,IIN,IOUT,ISAVE,J: Integer ' 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) FOR I = 2 TO N IF ABS(A1 - X(I)) < ACU THEN GOTO 10 A2 = X(I) IBAS2 = I Y2 = Y(I) GOTO 20 10 NEXT I 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 FOR I = 2 TO N DDD = Y(I) - (AAAA + BBBB * X(I)) 'D(I) = Sign(ONE, DDD) a = ONE: b = DDD: GOSUB 600: D(I) = Sign NEXT I TOT1 = ONE TOT2 = X(IBAS2) D(IBAS2) = -ONE FOR I = 2 TO N TOT1 = TOT1 + D(I) TOT2 = TOT2 + D(I) * X(I) NEXT I T = (A2 * TOT1 - TOT2) * DET IF ABS(T) < AONE THEN GOTO 50 DET = -DET GOTO 70 ' Main iterative loop begins 50 T = (TOT2 - A1 * TOT1) * DET IF ABS(T) < AONE THEN GOTO 160 IFLAG = 2 IOUT = IBAS2 X(IOUT) = A1 AAA = A1 BBB = A2 GOTO 80 60 T = (TOT2 - A2 * TOT1) * DET IF ABS(T) < AONE THEN GOTO 160 70 IFLAG = 1 BBB = A1 AAA = A2 IOUT = IBAS1 80 'RHO = Sign(ONE, T) a = ONE: b = T: GOSUB 600: RHO = Sign T = HALF * ABS(T) DET = DET * RHO ' Perform partial sort of ratios INEXT(IBAS1) = IBAS2 RATIO = BIG SUM = AHALF FOR I = 1 TO N DDD = (X(I) - AAA) * DET IF DDD * D(I) <= ACU THEN GOTO 120 TEST = (Y(I) - AAAA - BBBB * X(I)) / DDD IF TEST >= RATIO THEN GOTO 120 J = IBAS1 SUM = SUM + ABS(DDD) 90 ISAVE = ABS(INEXT(J)) IF TEST >= D(ISAVE) THEN GOTO 110 IF SUM < T THEN GOTO 100 SUBT = ABS((X(ISAVE) - AAA) * DET) IF SUM - SUBT < T THEN GOTO 100 SUM = SUM - SUBT 'D(ISAVE) = Sign(1, INEXT(J)) a = ONE: b = 1# * INEXT(J): GOSUB 600: D(ISAVE) = Sign INEXT(J) = INEXT(ISAVE) GOTO 90 100 J = ISAVE ISAVE = ABS(INEXT(J)) IF TEST < D(ISAVE) THEN GOTO 100 110 INEXT(I) = INEXT(J) 'INEXT(J) = ISign(I, Round(D^[I])) ia = I: ib = INT(D(I)): GOSUB 500: INEXT(J) = ISign D(I) = TEST IF SUM < T THEN GOTO 120 IIN = ABS(INEXT(IBAS1)) RATIO = D(IIN) 120 NEXT I ' Update basic indicators IIN = ABS(INEXT(IBAS1)) J = IIN 130 ISAVE = ABS(INEXT(J)) IF ISAVE = IBAS2 THEN GOTO 140 'ZZZ = ISign(1, INEXT(J)) ia = 1: ib = INEXT(J): GOSUB 500: ZZZ = ISign TOT1 = TOT1 - ZZZ - ZZZ TOT2 = TOT2 - TWO * ZZZ * X(ISAVE) D(ISAVE) = -ZZZ J = ISAVE GOTO 130 140 'ZZZ = ISign(1, INEXT(IBAS1)) ia = 1: ib = INEXT(IBAS1): GOSUB 500: ZZZ = ISign TOT1 = TOT1 - RHO - ZZZ TOT2 = TOT2 - RHO * BBB - ZZZ * X(IIN) D(IOUT) = -RHO ITER = ITER + 1 IF IFLAG = 1 THEN GOTO 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 GOTO 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 GOTO 50 ' Calculate optimal sum of absolute deviations 160 SAD = ZERO FOR I = 1 TO N D(I) = Y(I) - (AAAA + BBBB * X(I)) SAD = SAD + ABS(D(I)) NEXT I ALPHA = AAAA beta = BBBB RETURN 'end of file simlp.bas