! Last change: DOS 3 Aug 2000 6:06 pm ! *** copyright 2000 *** ! *** filename hlqest.f95 *** John F. Monahan ** ! ********************** PROGRAM PHLQEST IMPLICIT NONE REAL, DIMENSION(240) :: X REAL, DIMENSION(28920) :: Y REAL HLESTA,HLESTB REAL RAN,RANG,FSTMED,HLQEST ! 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,' H-L ESTIMATOR BY FSTMED=',F12.6) 23 FORMAT(' N=',I4,' H-L ESTIMATOR VIA HLQEST=',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='hlqest.out' ) ! initialize uniform generators HLESTA = RAN(5151917) HLESTB = 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 ! NN = 0 DO K = 1,N DO J = K,N NN = NN + 1 Y(NN) = ( X(K) + X(J) )/2. END DO ! LOOP ON J END DO ! LOOP ON K HLESTA = FSTMED(Y,NN) ! WRITE OUT H-L ESTIMATOR FOUND BY ENUMERATION WRITE(6,22) N,HLESTA ! ! COMPUTE H-L ESTIMATOR USING FUNCTION HLQEST ! HLESTB = HLQEST(X,N) WRITE(6,23) N,HLESTB IF( ABS( HLESTA - HLESTB ) .GT. 1.E-6 ) WRITE(6,29) ! ! now introduce some ties DO J = 1,N K = INT( 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 ! NN = 0 DO K = 1,N DO J = K,N NN = NN + 1 Y(NN) = ( X(K) + X(J) )/2. END DO ! LOOP ON J END DO ! LOOP ON K HLESTA = FSTMED(Y,NN) ! WRITE OUT H-L ESTIMATOR FOUND BY ENUMERATION WRITE(6,22) N,HLESTA ! ! COMPUTE H-L ESTIMATOR USING FUNCTION HLQEST ! HLESTB = HLQEST(X,N) WRITE(6,23) N,HLESTB IF( ABS( HLESTA - HLESTB ) .GT. 1.E-6 ) WRITE(6,29) ! done with this sample size END DO ! loop on i (N) STOP END PROGRAM PHLQEST REAL FUNCTION HLQEST(X,N) ! ! REAL FUNCTION HLQEST ! ! PURPOSE COMPUTES THE HODGES-LEHMANN LOCATION ESTIMATOR: ! MEDIAN OF ( X(I) + X(J) ) / 2 FOR 1 LE I LE J LE N ! ! USAGE RESULT = HLQEST(X,N) ! ! 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 * ! ! EXTERNAL ROUTINE ! RAN FUNCTION PROVIDING UNIFORM RANDOM VARIABLES ! IN THE INTERVAL (0,1) ! RAN REQUIRES A DUMMY INTEGER ARGUMENT ! ! NOTES HLQEST HAS AN EXPECTED TIME COMPLEXITY ON ! THE ORDER OF N * LG( N ) ! ! J F MONAHAN, APRIL 1982, DEPT OF STAT, N C S U, RALEIGH, N C 27650 ! FINAL VERSION JUNE 1983 ! RECODED MARCH, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN) :: X REAL AMN,AMX,AM,D INTEGER, DIMENSION(N) :: LB,RB,Q INTEGER SM,SQ,I,J,K,K1,K2,L,NN,LBI,RBI,MDLROW,IQ,METHOD REAL RAN ! ! TAKE CARE OF SPECIAL CASES: N=1 AND N=2 ! IF( N .EQ. 1 ) THEN HLQEST = X(1) RETURN END IF ! ( N .EQ. 1 ) IF( N .EQ. 2 ) THEN HLQEST = ( X(1) + X(2) )/2. 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 K1 = (NN+1)/2 K2 = (NN+2)/2 ! ! INITIALIZE LEFT AND RIGHT BOUNDS ! DO I = 1,N LB(I) = I RB(I) = N END DO ! LOOP ON I ! SM = NUMBER IN SET S AT STEP M SM = NN ! L = NUMBER OF PAIRS LESS THAN THOSE IN SET S AT STEP M L = 0 ! USE MEDIAN AS FIRST PARTITION ELEMENT METHOD = 1 ! *** MAIN LOOP BEGINNING *** DO ! WHAT TO USE AS PARTITION ELEMENT? SELECT CASE(METHOD) CASE(1) ! ! USE THE MEDIAN OF X(I)'S TO PARTITION ON THE FIRST STEP ! AM=X( (N+1)/2 ) + X( (N+2)/2 ) CASE(2) ! ! USE THE MIDRANGE OF SET S AS PARTITION ELEMENT WHEN TIES ARE LIKELY ! -- OR GET THE AVERAGE OF THE LAST 2 ELEMENTS ! AMX = X(1)+X(1) AMN = X(N)+X(N) DO I = 1,N ! SKIP THIS ROW IF NO ELEMENT IN IT IS IN SET S ON THIS STEP IF( LB(I) .LE. RB(I) ) THEN LBI = LB(I) ! GET THE SMALLEST IN THIS ROW AMN = MIN( AMN, X(LBI)+X(I) ) RBI=RB(I) ! GET THE LARGEST IN THIS ROW AMX = MAX( AMX, X(RBI)+X(I) ) END IF ! ( LB(I) .LE. RB(I) ) END DO ! LOOP ON I AM = (AMX+AMN)/2. ! BE CAREFUL TO CUT OFF SOMETHING -- ROUNDOFF CAN DO WIERD THINGS IF( (AM.LE.AMN) .OR. (AM.GT.AMX) ) AM=AMX ! UNLESS FINISHED, JUMP TO PARTITION STEP IF( (AMN.EQ.AMX) .OR. (SM.EQ.2) ) THEN ! ALL DONE IF ALL OF S IS THE SAME -OR- IF ONLY 2 ELEMENTS ARE LEFT HLQEST = AM/2. RETURN END IF ! (BOTH) CASE DEFAULT ! ! USE RANDOM ROW MEDIAN AS PARTITION ELEMENT K = INT(REAL(SM)*RAN(SM)) ! K IS A RANDOM INTEGER FROM O TO SM-1 DO I = 1,N J = I IF( K .LE. RB(I)-LB(I) ) EXIT K = K-RB(I)+LB(I)-1 END DO ! LOOP ON I ! J IS A RANDOM ROW --- NOW GET ITS MEDIAN MDLROW = (LB(J)+RB(J))/2 AM = X(J)+X(MDLROW) END SELECT ! (HOW TO PARTITION) ! ! ***** PARTITION STEP ***** ! ! USE AM TO PARTITION S0 INTO 2 GROUPS: THOSE .LT. AM, THOSE .GE. AM ! Q(I)= HOW MANY PAIRS (X(I)+X(J)) IN ROW I LESS THAN AM J = N ! START IN UPPER RIGHT CORNER SQ = 0 ! I COUNTS ROWS DO I = 1,N Q(I) = 0 ! HAVE WE HIT THE DIAGONAL ? IF( J .GE. I ) THEN ! SHALL WE MOVE LEFT ? DO WHILE ( J .GE. I ) D = X(I) + X(J) ! AVOID OPTIMIZATION/ARITHMETIC IF( D .LT. AM ) EXIT ! ANOMALIES -- BEWARE! J = J-1 END DO ! WHILE ( J .GE. I ) ! WE'RE DONE IN THIS ROW Q(I) = J-I+1 ! SQ = TOTAL NUMBER OF PAIRS LESS THAN AM SQ = SQ+Q(I) END IF ! ( J .GE. I ) END DO ! LOOP ON I ! ! *** FINISHED PARTITION --- START BRANCHING *** ! ! IF CONSECUTIVE PARTITIONS ARE THE SAME WE PROBABLY HAVE TIES IF( SQ .EQ. L ) THEN METHOD = 2 ! USE MIDRANGE ELSE ! ! ARE WE NEARLY DONE, WITH THE VALUES WE WANT ON THE BORDER? ! IF(WE NEED MAX OF THOSE .LT. AM -OR- MIN OF THOSE .GE. AM) EXIT ! IF( SQ .EQ. K2-1 ) EXIT IF( SQ .EQ. K1 ) EXIT ! ! THE SET S IS SPLIT, WHICH PIECE DO WE KEEP? ! SECOND= CUT OFF BOTTOM, EXIT = NEARLY DONE, FIRST= CUT OFF TOP ! IF(SQ .GT. K1) THEN ! ! NEW S = (OLD S) .INTERSECT. (THOSE .LT. AM) ! RESET RIGHT BOUNDS FOR EACH ROW DO I = 1,N RB(I) = I+Q(I)-1 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) END DO ! LOOP ON I END IF ! ( SQ .GT. K1 ) ! ! COUNT SM = NUMBER OF PAIRS STILL IN NEW SET S ! L = NUMBER OF PAIRS LESS THAN THOSE IN NEW SET S L = 0 SM = 0 DO I = 1,N L = L+LB(I)-I SM = SM+RB(I)-LB(I)+1 END DO ! LOOP ON I METHOD = 3 ! NORMAL ! END IF ! ( SQ .EQ. L ) IF( SM .EQ. 2 ) METHOD = 2 ! CAN ONLY GET TO 2 LEFT IF K1.NE.K2 -- GO GET THEIR AVERAGE END DO ! *** END OF MAIN LOOP *** ! ! FIND: MAX OF THOSE .LT. AM ! MIN OF THOSE .GE. AM ! AMN = X(N)+X(N) AMX = X(1)+X(1) DO I = 1,N IQ = Q(I) IF( IQ.GT.0) AMX = MAX( AMX, X(I)+X(I+IQ-1) ) IF( IQ.LT.N-I+1) AMN = MIN( AMN, X(I)+X(I+IQ) ) END DO ! LOOP ON I HLQEST = (AMN+AMX)/4. ! WE ARE DONE, BUT WHICH SITUATION ARE WE IN? IF( K1 .LT. K2 ) RETURN IF( SQ .EQ. K1 ) HLQEST = AMX/2. IF( SQ .EQ. K1-1 ) HLQEST = AMN/2. RETURN END FUNCTION HLQEST REAL FUNCTION FSTMED(X,N) ! FAST ALGORITHM FOR COMPUTING MEDIAN OF SAMPLE OF SIZE N ! ! ARGUMENTS ! X REAL VECTOR OF OBSERVATIONS ! N INTEGER NUMBER OF OBSERVATIONS ! ! *** REQUIRED SUBPROGRAMS *** ! SUBROUTINE PARTIT(X,LEFT,RGHT,NS,NT) PARTITIONS SET ! REAL FUNCTION RAN(IDUM) UNIFORM(0,1) RANDOM NUMBERS ! ! J F MONAHAN ! RECODED MARCH, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: X INTEGER J,K,KP,START,SPAN,NS,NT REAL RAN ! ! FIND MEDIAN -- K AND KP MAY OR MAY NOT BE DIFFERENT K = (N + 1) / 2 KP = (N + 2) / 2 ! INITIAL VALUES START = 1 SPAN = N ! RESTART HERE DO ! *** UNRESTRICTED DO *** ! RANDOM PARTITION ELEMENT (UNLESS TOO SHORT) NS = START IF( SPAN .GT. 2 ) NS = START + INT( RAN(SPAN)*REAL(SPAN) ) ! PARTITION X FROM START TO START+SPAN-1 WITH X(NS) CALL PARTIT(X,START,START+SPAN-1,NS,NT) ! NS = NUMBER SMALLER, NT = NUMBER TIED ! ! PARTITION MAKES TWO SUBLISTS ! ! ARE WE NEARLY DONE? IF( ( K .GE. START+NS-1 ) .AND. ( KP .LE. START+NS+NT ) ) EXIT ! LOOK AMONG SMALLER IF( K .LT. START + NS - 1 ) SPAN = NS ! LOOK AMONG LARGER IF( KP .GT. START + NS + NT ) THEN SPAN = SPAN - NS - NT START = START + NS + NT END IF ! ( KP .GT. START + NS + NT ) END DO ! *** END UNRESTRICTED DO *** ! NEARLY DONE -- MAY HAVE TO FIND MIN OR MAX IF( KP .EQ. START+NS+NT ) THEN ! FIND MIN OF LARGERS FSTMED = X(START+NS+NT) SPAN = SPAN - NS - NT DO J = 1,SPAN FSTMED = MIN( FSTMED, X(START+NS+NT+J-1) ) END DO ! LOOP ON J IF( K .NE. KP ) FSTMED = ( FSTMED + X(START+NS) ) / 2. RETURN END IF ! ( KP .EQ. START+NS+NT ) ! IF( K .EQ. START+NS-1 ) THEN ! FIND MAX OF SMALLERS FSTMED = X(START) DO J = 1,NS FSTMED = MAX( FSTMED, X(START+J-1) ) END DO ! LOOP ON J IF( K .NE. KP ) FSTMED = ( FSTMED + X(START+NS) ) / 2. RETURN END IF ! ( K .EQ. START+NS-1 ) ! DONE -- MEDIAN IN MIDDLE FSTMED = X(START+NS) RETURN END FUNCTION FSTMED SUBROUTINE PARTIT(X,LEFT,RGHT,NS,NT) ! ALGORITHM FOR PARTITIONING LIST OF KEYS K FROM X(LEFT) TO X(RGHT) ! INTO THREE GROUPS: THOSE < X*, THOSE = X*, THOSE > X* ! ! ARGUMENTS ! X IN REAL LIST OF KEYS TO BE PARTITIONED ! OUT PARTITIONED LIST ! X(LEFT) ... X(LEFT+NS-1) < X* ! X(LEFT+NS) ... X(LEFT+NS+NT-1) = X* ! X(LEFT+NS+NT) ... X(RGHT) > X* ! ! LEFT IN INTEGER LEFT ENDPOINT OF LIST ! RGHT IN INTEGER RIGHT ENDPOINT OF LIST ! NS IN INTEGER POINTS TO PARTITION ELEMENT X* = X(NS) (IN) ! OUT NUMBER IN LIST < X* ! NT OUT INTEGER NUMBER IN LIST = X* ! IMPLICIT NONE INTEGER, INTENT(IN) :: LEFT,RGHT INTEGER, INTENT(IN OUT) :: NS INTEGER, INTENT(OUT) :: NT REAL, DIMENSION(RGHT) :: X REAL XCUT,XX INTEGER I,J ! USE X(NS) TO SPLIT X XCUT = X(NS) ! DO NOTHING IF THERE'S NOTHING TO DO NS = 0 NT = 1 IF( RGHT - LEFT .LE. 0 ) RETURN ! FIRST PARITITION TO TWO GROUPS: LT AND GE ! ! INITIALIZE POINTERS I (LEFT) AND J (RIGHT) I = LEFT J = RGHT DO WHILE( I .LE. J ) DO WHILE( X(I) .LT. XCUT ) I = I + 1 IF( I .GT. J ) EXIT ! LAST ONE LT IS I-1 END DO ! WHILE ( X(I) .LT. XCUT ) IF( I .GT. J ) EXIT ! FINISH ! HAVE X(I) OUT OF PLACE DO WHILE ( X(J) .GE. XCUT ) J = J - 1 IF( I .GE. J ) EXIT END DO ! WHILE ( X(J) .GE. XCUT ) IF( I .GE. J ) EXIT ! HAVE X(J) OUT OF PLACE ! EXCHANGE I AND J XX = X(I) X(I) = X(J) X(J) = XX I = I + 1 J = J - 1 END DO ! FIRST MAIN WHILE ! HALFWAY DONE -- X(I-1) IS END OF SMALLER NS = I - LEFT ! NOW SEPARATE THE TIED FROM THE LARGER ! ! ! INITIALIZE POINTERS I (LEFT) AND J (RIGHT) J = RGHT DO WHILE( I .LE. J ) DO WHILE( X(I) .EQ. XCUT ) I = I + 1 IF( I .GT. J ) EXIT ! LAST ONE EQ IS I-1 END DO ! WHILE ( X(I) .EQ. XCUT ) IF( I .GT. J ) EXIT ! FINISH ! HAVE X(I) OUT OF PLACE DO WHILE( X(J) .GT. XCUT ) J = J - 1 IF( I .GE. J ) EXIT END DO ! WHILE ( X(J) .GT. XCUT ) IF( I .GE. J ) EXIT ! HAVE X(J) OUT OF PLACE ! EXCHANGE I AND J XX = X(I) X(I) = X(J) X(J) = XX I = I + 1 J = J - 1 END DO ! SECOND MAIN WHILE NT = I - LEFT - NS RETURN END SUBROUTINE PARTIT 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 hlqest.f95 *********************