! Last change: DOS 3 Aug 2000 5:46 pm ! *** copyright 2000 *** ! *** filename abdftr.f95 *** John F. Monahan ** ! ********************** PROGRAM PABDFTR IMPLICIT NONE REAL, DIMENSION(240) :: X REAL, DIMENSION(28920) :: Y REAL SESTA,SESTB REAL RAN,RANG,ABDFTR ! RANG JUST A SECOND COPY INTEGER I,J,K,N,NN ! 20 FORMAT(//' Test Problems with Sample Size',i8) 21 FORMAT(8F9.4) 22 FORMAT(' N=',I4,' SCALE ESTIMATOR BY SORTING=',F12.6) 23 FORMAT(' N=',I4,' SCALE ESTIMATOR VIA ABDFTR=',F12.6) 25 FORMAT(/' Simple Case -- Mixed Signs, N=',i8) 27 FORMAT(/' Mixed Signs and Some Ties, N=',i8) 29 FORMAT(' ESTIMATORS ARE DIFFERENT' ) ! open( unit=6, file='abdftr.out' ) ! initialize uniform generators SESTA = RAN(5151917) SESTB = RANG(5151917) ! let's check out these twenty cases DO I = 1,20 N = I IF( I .GT. 10 ) N = I*(I-8) WRITE(6,20) N ! now mixed positive and negative and smaller DO J = 1,N X(J) = (RANG(J)-.3)/100. END DO ! loop on j WRITE(6,25) N ! ! SORT THE X(I)'S CALL HSORT(X,N) IF( N .LE. 10 ) WRITE(6,21) (X(J),J=1,N) ! ! STRAIGHTFORWARD METHOD: GET ALL N(N-1)/2 PAIRS AND SORT ! NN = 0 DO J = 2,N DO K = 1,J-1 NN = NN + 1 Y(NN) = ( X(J) - X(K) ) END DO ! LOOP ON J END DO ! LOOP ON K CALL HSORT(Y,NN) K = INT( REAL(NN)*RANG(I) ) + 1 ! CHOOSE K SESTA = 0. DO J = 1,K SESTA = SESTA + Y(J) END DO ! LOOP ON J ! WRITE OUT SCALE ESTIMATOR FOUND BY ENUMERATION WRITE(6,22) N,SESTA ! ! COMPUTE SCALE ESTIMATOR USING FUNCTION ABDFTR ! SESTB = ABDFTR(X,N,K) WRITE(6,23) N,SESTB IF( ABS( SESTA - SESTB ) .GT. 1.E-5*SESTA ) WRITE(6,29) ! ! now introduce some ties DO J = 1,N K = 8.*(RANG(J)-.3) X(J) = REAL(K)/10. END DO ! loop on j WRITE(6,27) N ! ! SORT THE X(I)'S CALL HSORT(X,N) IF( N .LE. 10 ) WRITE(6,21) (X(J),J=1,N) ! ! STRAIGHTFORWARD METHOD: GET ALL N(N-1)/2 PAIRS AND SORT ! NN = 0 DO J = 2,N DO K = 1,J-1 NN = NN + 1 Y(NN) = ( X(J) - X(K) ) END DO ! LOOP ON J END DO ! LOOP ON K CALL HSORT(Y,NN) K = INT( REAL(NN)*RANG(I) ) + 1 ! CHOOSE K SESTA = 0. DO J = 1,K SESTA = SESTA + Y(J) END DO ! LOOP ON J ! WRITE OUT SCALE ESTIMATOR FOUND BY ENUMERATION WRITE(6,22) N,SESTA ! ! COMPUTE SCALE ESTIMATOR USING FUNCTION ABDFTR ! SESTB = ABDFTR(X,N,K) WRITE(6,23) N,SESTB IF( ABS( SESTA - SESTB ) .GT. 1.E-5*SESTA ) WRITE(6,29) ! done with this sample size end do ! loop on i (N) STOP END PROGRAM PABDFTR REAL FUNCTION ABDFTR(X,N,K) ! ! REAL FUNCTION ABDFTR ! ! PURPOSE COMPUTES THE SCALE ESTIMATOR: SUM OF THE K SMALLEST ! VALUES OF ( X(J) - X(I) ) / 2 FOR 1 LE I LT J LE N ! ! USAGE RESULT = HLQEST(X,N,K) ! ! ARGUMENTS X REAL ARRAY OF OBSERVATIONS (INPUT) ! * VALUES OF X MUST BE IN NONDECREASING ORDER * ! ! N INTEGER NUMBER OF OBSERVATIONS (INPUT) ! * N MUST NOT BE LESS THAN 1 * ! ! K NUMBER OF PAIRS OF DIFFERENCES TO BE ADDED ! ! EXTERNAL ROUTINE ! RAN FUNCTION PROVIDING UNIFORM RANDOM VARIABLES ! IN THE INTERVAL (0,1) ! RAN REQUIRES A DUMMY INTEGER ARGUMENT ! ! NOTES ABDFTR HAS AN EXPECTED TIME COMPLEXITY ON ! THE ORDER OF N * LG( N ) ! ! J F MONAHAN, AUGUST 1986, DEPT OF STAT, N C S U, RALEIGH, N C 27650 ! MODIFICATION OF ALGORITHM 616:HLQEST, ACM TOMS (1984) V 10, P265-70 ! RECODED MARCH, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N,K REAL, DIMENSION(N), INTENT(IN) :: X REAL AMN,AMX,AM,D INTEGER, DIMENSION(N) :: LB,RB,Q INTEGER SM,SQ,I,J,L,NN REAL RAN ! ! TAKE CARE OF SPECIAL CASES: N=1 AND N=2 ! IF( N .EQ. 1 ) THEN ABDFTR = 0. RETURN END IF ! ( N .EQ. 1 ) IF( N .EQ. 2 ) THEN ABDFTR = X(2) - X(1) RETURN END IF ! ( N .EQ. 2 ) ! ! FIND THE TOTAL NUMBER OF PAIRS (NN) AND THE MEDIAN(S) (K1,K2) NEEDED ! NN = (N*(N-1))/2 ! ! INITIALIZE LEFT AND RIGHT BOUNDS ! DO I = 1,N LB(I) = I + 1 RB(I) = N END DO ! LOOP ON I ! SM = NUMBER IN SET S AT STEP M SM = NN ! *** MAIN LOOP BEGINNING *** DO ! ! USE RANDOM ELEMENT AS PARTITION ELEMENT L = INT(REAL(SM)*RAN(SM)) ! K IS A RANDOM INTEGER FROM O TO SM-1 DO I = 1,N J = I IF( L .LE. RB(I)-LB(I) ) EXIT L = L-RB(I)+LB(I)-1 END DO ! LOOP ON I ! J IS A RANDOM ROW --- NOW GET ITS MEDIAN L = LB(J) + L AM = X(L) - X(J) ! ! ***** PARTITION STEP ***** ! ! USE AM TO PARTITION S0 INTO 2 GROUPS: THOSE .LT. AM, THOSE .GE. AM ! Q(I)= HOW MANY PAIRS (X(J)-X(I)) IN ROW I LESS THAN AM J = 1 ! START IN UPPER LEFT CORNER SQ = 0 ! I COUNTS ROWS DO I = 1,N J = MAX( J, I+1 ) ! HAVE WE HIT THE DIAGONAL ? IF( J .LE. N ) THEN ! SHALL WE MOVE RIGHT ? DO WHILE ( J .LE. N ) D = X(J) - X(I) ! AVOID OPTIMIZATION/ARITHMETIC IF( D .GE. AM ) EXIT ! ANOMALIES -- BEWARE ! J = J + 1 END DO ! WHILE ( J .LE. N ) END IF ! ( J .LE. N ) ! WE'RE DONE IN THIS ROW Q(I) = J - I - 1 ! SQ = TOTAL NUMBER OF PAIRS LESS THAN AM SQ = SQ + Q(I) END DO ! LOOP ON I ! ! *** FINISHED PARTITION --- START BRANCHING *** ! ! ! THE SET S IS SPLIT, WHICH PIECE DO WE KEEP? ! SECOND= CUT OFF BOTTOM, FIRST= CUT OFF TOP ! IF(SQ .GT. K) THEN ! ! NEW S = (OLD S) .INTERSECT. (THOSE .LT. AM) ! RESET RIGHT BOUNDS FOR EACH ROW DO I = 1,N RB(I) = I+Q(I) Q(I) = 0 END DO ! LOOP ON I ELSE ! NEW S = (OLD S) .INTERSECT. (THOSE .GE. AM) ! RESET LEFT BOUNDS FOR EACH ROW DO I = 1,N LB(I) = I+Q(I)+1 Q(I) = 0 END DO ! LOOP ON I END IF ! ( SQ .GT. K ) ! ! COUNT SM = NUMBER OF PAIRS STILL IN NEW SET S ! AMX = 0. AMN = X(N) - X(1) SM = 0 DO I = 1,N SM = SM+RB(I)-LB(I)+1 IF( RB(I) .GE. LB(I) ) THEN J = RB(I) AMX = MAX(AMX, X(J) - X(I) ) J = LB(I) AMN = MIN(AMN, X(J) - X(I) ) END IF ! ( RB(I) .GE. LB(I) ) END DO ! LOOP ON I ! IF( SM .EQ. 1 ) EXIT IF( AMN .EQ. AMX ) EXIT ! END DO ! *** END OF MAIN LOOP *** ! ! FINISH UP BY ADDING UP EVERYTHING ! ESTIMATOR CAN BE WRITTEN AS WEIGHTED ! SUM OF OBSERVATIONS, GET Q(I) WEIGHTS ! ! FIRST GET THE WEIGHTS FOR THE OBSERVATIONS L = 0 J = 1 DO I = 1,N ! FIRST SUBTRACT X(I) FOR EVERY ELEMENT IN ROW I Q(I) = Q(I) + (I+1) - LB(I) L = L + LB(I) - (I+1) ! NOW ADD X(J) FOR EVERY ELEMENT IN COLUMN J DO WHILE ( J .LT. LB(I) ) Q(J) = Q(J) + J - I J = J + 1 END DO ! WHILE END DO ! LOOP ON I ! ADD THEM UP TO GET ESTIMATOR ABDFTR = 0. DO I = 1,N ABDFTR = ABDFTR + Q(I)*X(I) END DO ! LOOP ON I ! DON'T FORGET THE TIES IF( K .NE. L ) ABDFTR = ABDFTR + (K-L)*AMN RETURN END FUNCTION ABDFTR REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN REAL FUNCTION RANG(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RANG = REAL(IX)*RP RETURN END FUNCTION RANG SUBROUTINE HSORT(K,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: K REAL KK INTEGER I,L,NCUR ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HSORT ! *** end of filename abdftr.f95 *********************