PROCEDURE ROBSCALE MYSERIES NBEG NEND
*
* ROBSCALE.SRC COMPUTES A ROBUST ALTERNATIVE TO THE STANDARD
* DEVIATION AS A MEASURE OF SCALE (DISPERSION). THE STANDARD
* DEVIATION IS EASILY DISTORTED BY OUTLIERS AND MAY THEREFORE
* BE QUITE MISLEADING. ROBSCALE.SRC COMPUTES THE Qn STATISTIC
* OF ROUSSEEUUW AND CROUX, WHICH IS VERY ROBUST AND ALSO HAS
* GOOD EFFICIENCY (RELATIVE TO THE STANDARD DEVIATION) WHEN THE
* DATA ARE IN FACT UNCONTAMINATED GAUSSIAN RANDOM VARIABLES.
* (P. J. ROUSSEEUW AND C. CROUX, "ALTERNATIVES TO THE MEDIAN
* ABSOLUTE DEVIATION," JOURNAL OF THE AMERICAN STATISTICAL
* ASSOCIATION, DEC 1993, VOL. 88, NO. 424, PP. 1273-1283).
* SPECIFICALLY, IF THE SERIES X CONTAINS N OBSERVATIONS, THEN
* THERE ARE N(N-1)/2 DIFFERENCES OF THE FORM ABS( X(I)-X(J) ),
* I < J. WHEN THESE DIFFERENCES ARE RANKED, Qn IS APPROXIMATELY
* THE 25TH PERCENTILE. THE COST OF COMPUTING Qn GROWS RAPIDLY
* AS N INCREASES, BUT THE ROUSSEEUW-CROUX ALGORITHM IS VERY
* FAST AND NEEDS RELATIVELY LITTLE COMPUTER MEMORY AND DISK
* SPACE.
*
* MYSERIES IS THE UNIVARIATE DATA FROM WHICH Qn IS TO BE
* COMPUTED. THE FIRST OBSERVATION IS NBEG AND THE LAST
* OBSERVATION IS NEND.
*
TYPE SERIES MYSERIES
TYPE INTEGER NBEG NEND
LOCAL INTEGER I J JJ N H KK KNEW JHELP NL NR SUMQ SUMP FOUND
LOCAL INTEGER WTOTAL WREST WLEFT WMID WRIGHT NN KCAND
LOCAL VEC WORK ACAND Y
LOCAL VEC[INTEGER] LEFT RIGHT WEIGHT
LOCAL VEC[INTEGER]  Q P IWCAND
LOCAL SERIES XSERIES WORK1 WORK2
COMPUTE N = NEND-NBEG+1
DIM WORK(N) ACAND(N) Y(N)
DIM LEFT(N) RIGHT(N) WEIGHT(N)
DIM Q(N) P(N) IWCAND(N)
COMPUTE H=N/2+1
COMPUTE KK=H*(H-1)/2
SET XSERIES NBEG NEND = MYSERIES
ORDER XSERIES NBEG NEND
COMPUTE J = 0
DO I = NBEG,NEND
COMPUTE J=J+1
COMPUTE Y(J) = XSERIES(I)
END DO I
DO I = 1,N
COMPUTE LEFT(I)=N-I+2
COMPUTE RIGHT(I)=N
END DO I
COMPUTE JHELP=N*(N+1)/2
COMPUTE KNEW=KK+JHELP
COMPUTE NL=JHELP
COMPUTE NR=N*N
COMPUTE FOUND = 0
DO MGM = 1,1
:200 IF((NR-NL.LE.N) .OR. (FOUND.EQ.1))
BRANCH 1010
COMPUTE J=1
DO I =2,N
IF (LEFT(I).LE.RIGHT(I))
   {
   COMPUTE WEIGHT(J)=RIGHT(I)-LEFT(I)+1
   COMPUTE JHELP=LEFT(I)+WEIGHT(J)/2
   COMPUTE WORK(J)=Y(I)-Y(N+1-JHELP)
   COMPUTE J=J+1
   }
END DO I
COMPUTE NN=J-1
COMPUTE WTOTAL=0
DO I=1,NN
COMPUTE WTOTAL=WTOTAL+WEIGHT(I)
END DO I
COMPUTE WREST=0
:100 SET WORK1 1 NN = WORK(T)
ORDER WORK1 1 NN
COMPUTE TRY = WORK1(NN/2+1)
COMPUTE WLEFT=0
COMPUTE WMID=0
COMPUTE WRIGHT=0
DO I = 1,NN
IF (WORK(I) .LT. TRY)
COMPUTE WLEFT=WLEFT+WEIGHT(I)
IF (WORK(I) .GT. TRY)
COMPUTE WRIGHT=WRIGHT+WEIGHT(I)
IF (WORK(I) .EQ. TRY)
COMPUTE WMID=WMID+WEIGHT(I)
END DO I
IF ((2*WREST+2*WLEFT).LE.WTOTAL)
BRANCH 5050
COMPUTE KCAND=0
DO I = 1,NN
IF (WORK(I).LT.TRY)
   {
   COMPUTE KCAND=KCAND+1
   COMPUTE ACAND(KCAND)=WORK(I)
   COMPUTE IWCAND(KCAND)=WEIGHT(I)
   }
END DO I
COMPUTE NN=KCAND
BRANCH 6060
:5050 IF ((2*WREST+2*WLEFT+2*WMID).GT.WTOTAL)
   {
   COMPUTE WHIMED=TRY
   BRANCH 900
   }
COMPUTE KCAND=0
DO I=1,NN
IF (WORK(I).GT.TRY)
   {
   COMPUTE KCAND=KCAND+1
   COMPUTE ACAND(KCAND)=WORK(I)
   COMPUTE IWCAND(KCAND)=WEIGHT(I)
   }
END DO I
COMPUTE NN=KCAND
COMPUTE WREST=WREST+WLEFT+WMID
:6060 DO I = 1,NN
COMPUTE WORK(I)=ACAND(I)
COMPUTE WEIGHT(I) = IWCAND(I)
END DO I
BRANCH 100
:900 COMPUTE TRIAL= WHIMED
COMPUTE J=0
DO I=N,1,-1
:45 IF ((J.LT.N) .AND. ((Y(I)-Y(N-J)) .LT. TRIAL))
   {
   COMPUTE J=J+1
   BRANCH 45
   }
COMPUTE P(I)=J
END DO I
COMPUTE J=N+1
DO I =1,N
:55 IF((Y(I)-Y(N-J+2)).GT.TRIAL)
  {
  COMPUTE J=J-1
  BRANCH 55
  }
COMPUTE Q(I)=J
END DO I
COMPUTE SUMP=0
COMPUTE SUMQ=0
DO I=1,N
COMPUTE SUMP=SUMP+P(I)
COMPUTE SUMQ=SUMQ+Q(I)-1
END DO I
IF (KNEW.GT.SUMP)
BRANCH 2020
DO I=1,N
COMPUTE RIGHT(I)=P(I)
END DO I
COMPUTE NR=SUMP
BRANCH 3030
:2020 IF (KNEW.LE.SUMQ)
BRANCH 4040
DO I=1,N
COMPUTE LEFT(I)=Q(I)
END DO I
COMPUTE NL=SUMQ
BRANCH 3030
:4040 COMPUTE QN=TRIAL
COMPUTE FOUND=1
:3030 BRANCH 200
:1010 IF FOUND .EQ. 0
   {
   COMPUTE J=1
   DO I =2,N
   IF (LEFT(I) .LE. RIGHT(I))
           {
           DO JJ=LEFT(I),RIGHT(I)
           COMPUTE WORK(J)=Y(I)-Y(N-JJ+1)
           COMPUTE J=J+1
           END DO JJ
           }
     END DO I
     SET WORK2 1 J-1 = WORK(T)
     ORDER WORK2 1 J-1
     COMPUTE QN = WORK2(KNEW-NL)
     }
IF (N .GT. 9)
BRANCH 7070
IF (N.EQ.2)
COMPUTE DN=0.399
IF (N.EQ.3)
COMPUTE DN=0.994
IF (N.EQ.4)
COMPUTE DN=0.512
IF (N.EQ.5)
COMPUTE DN=0.844
IF (N.EQ.6)
COMPUTE DN=0.611
IF (N.EQ.7)
COMPUTE DN=0.857
IF (N.EQ.8)
COMPUTE DN=0.669
IF (N.EQ.9)
COMPUTE DN=0.872
BRANCH 8080
:7070 COMPUTE RN=N
COMPUTE HARPO = %FRAC(RN/2.0)
IF HARPO .NE. 0.000
COMPUTE DN=N/(N+1.4)
IF HARPO .EQ. 0.000
COMPUTE DN=N/(N+3.8)
:8080 COMPUTE QN = DN*2.2219*QN
END MGM
END

