! Last change: DOS 3 Aug 2000 6:36 pm ! *** copyright 2000 *** ! *** filename wlcx1s.f95 *** John F. Monahan ** ! ********************** program wlcx1s ! DEMONSTRATION OF WLCX1S ! COMPARISON OF PAGANO-TRITCHLER METHOD FOR COMPUTING DIST OF RANK TESTS ! AND WLCXPT -- WHICH COMPUTES THE DISTRIBUTION ! OF THE ONE-SAMPLE WILCOXON STATISTIC WHEN TIES ARE PRESENT ! implicit none complex, dimension(256) :: phi real t,s real, parameter :: twopi = 6.2831853 integer, dimension(821) :: c integer, dimension(20) :: ties integer j,k,length ! ! INTERFACE BLOCK INTERFACE SUBROUTINE WLCXPT(T,Q,C,LNEW) INTEGER, INTENT(IN) :: Q INTEGER, DIMENSION(:), INTENT(IN) :: T INTEGER, DIMENSION(:), INTENT(OUT) :: C INTEGER, INTENT(OUT) :: LNEW END SUBROUTINE WLCXPT END INTERFACE ! 20 FORMAT(' Sample Size=',i3,4x,'Max Value of Statistic',I4 & & //' Probability Distribution * 2**N') 21 FORMAT(2X,25I3) 22 format(i4,3f12.6,4x,f8.6) 23 format(/9x,'Prob from',3x,'Real Part',3x,'Imag Part',4x,'Sum of',& & /' I',6x,'WLCXPT',6x,'of Phi',6x,'of Phi',6x,'Probs') ! data ties/20*1/ ! open( unit=6, file='wlcx1s.out' ) ! ! DO SAMPLE SIZE 7 WITHOUT TIES TO CHECK CALL WLCXPT(TIES,7,C,LENGTH) WRITE(6,20) 7,LENGTH WRITE(6,21) (C(J),J=1,LENGTH) ! now do it according to Pagano and Trichler ! do k = 1,64 t = -twopi*real(k-1)/64. ! get characteristic function at this t ! add exp(i*sum) to get usual definition of statistic phi(k) = cmplx( cos(28.*t), sin(28.*t) ) do j = 1,7 phi(k) = phi(k) * cmplx( cos(t*float(j)) , 0. ) end do ! loop on j end do ! loop on k ! invert characteristic function call fft2n(phi,6,-1.) ! write out results write(6,23) s = 0. do k = 1,length ! probability computed using WLCXPT t = c(k)/real(2**7) ! s gives cumulative prob s = s + real( phi(k) ) write(6,22) k,t,phi(k),s end do ! loop on k stop end program wlcx1s SUBROUTINE WLCXPT(T,Q,C,LNEW) ! COMPUTES THE NULL DISTRIBUTION OF THE ONE-SAMPLE WILCOXON STATISTI! ! WHEN TIES MAY BE PRESENT ! ! T INTEGER VECTOR INPUT HOLDS PATTERN OF TIES (NO TIES T(I) = 1 ) ! * NOTE * NO MORE THAN 50 TIES AT ONE LEVEL ! ! Q INTEGER INPUT NUMBER OF TIED GROUPS (NO TIES Q = NOBS) ! ! C INTEGER VECTOR OUTPUT HOLDS PROBABILITY VALUES OF DISTRIBUTION ! DISPLACED BY 1 ( PR(W=0) = C(1)/2**N ) ! * NOTE * NORMALIZATION CONSTANT IS 2**N ! ! LNEW INTEGER OUTPUT GIVES LENGTH OF DISTRIBUTION VECTOR C ! * NOTE * BECAUSE OF TIES, INCREMENTS ARE ! BY HALF INTEGER, SO MAX VALUE IS (LNEW-1)/2 ! ! J F MONAHAN (AUGUST,1989) DEPT OF STATISTICS, NCSU, RALEIGH, NC 27695 ! RECODED MARCH 2000 FOR FORTRAN 95 ! INTEGER, INTENT(IN) :: Q INTEGER, DIMENSION(:), INTENT(IN) :: T INTEGER, DIMENSION(:), INTENT(OUT) :: C INTEGER, INTENT(OUT) :: LNEW INTEGER, DIMENSION(51) :: BC INTEGER LNOLD,I,J,K,TJ,ST,KD,STEP,S,TJP1 ! ! INITIALIZE PROBABILITY VECTOR TO ZEROS LNEW = 1 ST = 0 DO J = 1,Q LNEW = LNEW + T(J)*(2*ST + T(J) + 1) ST = ST + T(J) END DO ! LOOP ON J DO K = 1,LNEW C(K) = 0 END DO ! LOOP ON K ! DO FIRST GROUP AND DO TIES AS BINOMIALS TJ = T(1) BC(1) = 1 DO I = 1,TJ BC(I+1) = ( (TJ-I+1)*BC(I) )/I END DO ! LOOP ON I C(1) = 1 STEP = TJ+1 DO I = 1,TJ KD = I*STEP + 1 C(KD) = BC(I+1) END DO ! LOOP ON I LNOLD = KD ! LNOLD = CURRENT LENGTH, ST = RUNNING SUM OF TIES ST = T(1) IF( Q .EQ. 1 ) RETURN DO J =2,Q TJ = T(J) TJP1 = TJ + 1 ! TABLE BINOMIAL COEFFICIENTS BC(1) = 1 DO I = 1,TJ BC(I+1) = ( (TJ-I+1)*BC(I) )/I END DO ! LOOP ON I STEP = 2*ST + TJ + 1 LNEW = LNOLD + TJ * STEP ! STEP THROUGH EACH POSSIBLE ONE FROM BIGGEST DOWN DO K = 1,LNEW S = 0. DO I = 1,TJP1 ! CONTRIBUTORS TO LNEW-K+1 ARE DISPLACED BY STEP UNITS KD = LNEW - K + 1 - (I-1)*STEP IF ( (KD .GE. 1 ) .AND. (KD .LE.LNOLD) ) S = S + C(KD)*BC(I) END DO ! LOOP ON I C(LNEW-K+1) = S END DO ! LOOP ON K ST = ST + TJ LNOLD = LNEW END DO ! LOOP ON J RETURN END SUBROUTINE WLCXPT SUBROUTINE FFT2N(A,N,SGN) ! ! SUBROUTINE FFT2N -- DISCRETE FOURIER TRANSFORM ! ! A COMPLEX VECTOR INPUT: VECTOR TO BE TRANSFORMED ! OUTPUT: TRANSFORMED VECTOR ! N INTEGER INPUT: LOG(BASE 2) OF LENGTH OF VECTOR A ! LENGTH OF A IS 2 ** N ! SGN REAL INPUT: INDICATES FORWARD OR INVERSE TRANSFO ! SGN = 1. THEN FORWARD TRANSFORM ! SGN = -1. THEN INVERSE TRANSFORM ! ! DISCRETE FOURIER TRANSFORM OF A VECTOR A(J) IS GIVEN BY ! NN-1 ! T( A ) = SUM A(J) EXP( - I SGN 2 PI J K / NN ) ! J K=0 ! ! WHERE I = SQRT( -1 ) ! NN = 2 ** N ! PI = 3.1415923565... ! ! NOTE THAT THE INDEXING OF A IS SHIFTED BY ONE -- ! A(0) STORED IN A(1), ... , A( NN-1 ) STORED IN A( NN ) ! ! J F MONAHAN (NOV,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695-8203 ! RECODED DECEMBER 1999, APRIL 2000 FOR FORTRAN 95 ! ! THIS IS A 'POWER OF 2' ALGORITHM BASED ON THE PRESENTATION IN ! ! A. V. AHO, J. E. HOPCROFT, AND J. D. ULLMAN (1974) THE DESIGN AND ! ANALYSIS OF COMPUTER ALGORITHMS, ADDISON-WESLEY, READING. ! ! ! BEGIN WITH DECLARATIONS IMPLICIT NONE INTEGER, INTENT(IN) :: N COMPLEX, DIMENSION(1:N), INTENT(IN OUT) :: A REAL, INTENT(IN) :: SGN REAL C,S COMPLEX, DIMENSION(25) :: WK COMPLEX WS,ALEFT,ARIGHT LOGICAL CARRY,NEW LOGICAL, DIMENSION(25) :: RBITS INTEGER I,J,L,LM,LL,M,MM,M2,NN ! ! COMPUTE WK(I) = EXP( - 2 PI IMAG SGN / 2**I ) ! ! SPECIAL VALUES FIRST WK(1) = CMPLX( -1., 0. ) WK(2) = CMPLX( 0., -SGN ) ! SKIP IF N IS ONLY 1 OR 2 IF( N .GT. 2 ) THEN ! AVOID SINES AND COSINES -- JUST USE SQUARE ROOTS C = 0. S = 1. DO I = 3,N C = SQRT( ( 1. + C ) / 2. ) S = S / ( 2.*C ) WK(I) = CMPLX( C, -SGN*S ) END DO ! LOOP ON I END IF ! ( N .GT. 2 ) ! ! THIS IS THE MAIN LOOP OF THE ALGORITHM ! DO MM = 1,N ! CONVERT TO REVERSE LOOP M = N - MM ! M STEPS DOWN FROM N-1 TO 0 ! ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I ! M2 M2 = 2**M ! LM CONTROLS THE NUMBER OF SUBLOOPS LM = 2**(MM-1) ! DO LL = 1,LM ! L IS THE MAJOR INDEX VARIABLE L = 2*M2*(LL-1) ! INITIALIZE TO GET W ** REV( L / M2 ) WS = CMPLX(1., 0.) CARRY = .TRUE. ! START AT 2 (SHIFT) SINCE L/M2 IS ALWAYS EVEN DO I = 2,N NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) WS = WS * WK(I) END DO ! LOOP ON I ! WS IS NOW W ** REV( L / M2 ) ! NOW DO THE REAL COMPUTATION DO J = 1,M2 ALEFT = A( L + J ) ARIGHT = WS*A( L + J + M2 ) A( L + J ) = ALEFT + ARIGHT A( L + J + M2 ) = ALEFT - ARIGHT END DO ! LOOP ON J ! DONE WITH LOOP ON L END DO ! LOOP ON LL ! DONE WITH MAIN LOOP ON M END DO ! LOOP ON MM ! ! SCRAMBLE THE NUMBERS -- SWITCH A(J) AND A( REV(J) ) ! BUT DON'T DO IT TWICE ! INITIALIZE DO I = 1,N RBITS(I) = .TRUE. END DO ! LOOP ON I NN = 2 ** N ! LOOP THROUGH ALL THE NUMBERS DO LL = 1,NN ! DO THE BIT REVERSAL LM = 0 M2 = NN CARRY = .TRUE. ! START AT 1 HERE -- USE ALL THE NUMBERS DO I = 1,N M2 = M2 / 2 NEW = CARRY .OR. RBITS(I) CARRY = CARRY .AND. RBITS(I) RBITS(I) = NEW .AND. ( .NOT. CARRY ) IF( RBITS(I) ) LM = LM + M2 END DO ! LOOP ON I ! DON'T SWITCH TWICE IF( LL-1 .LT. LM ) THEN ALEFT = A(LL) ARIGHT = A(LM+1) A(LL) = ARIGHT A(LM+1) = ALEFT END IF ! ( LL-1 .LT. LM ) END DO ! LOOP ON LL ! IF INVERSE TRANSFORM, REMEMBER TO DIVIDE BY NN IF( SGN .GT. 0. ) RETURN C = REAL(NN) DO I = 1,NN A(I) = A(I) / C END DO ! LOOP ON I RETURN END SUBROUTINE FFT2N ! *** end of filename wlcx1s.f95 *********************