! Last change: DOS 3 Aug 2000 6:26 pm ! *** copyright 2000 *** ! *** filename t2perm.f95 *** John F. Monahan ** ! ********************** PROGRAM T2PERM ! PERMUTATION TWO SAMPLE T-TEST -- GET THE (ONE-SIDED) P-VALUE ! IMPLICIT NONE REAL, DIMENSION(20) :: X,Y REAL, DIMENSION(40) :: BOTH REAL TORIG,T REAL RBICOF INTEGER, DIMENSION(20) :: LIST INTEGER M,N,I,KL,KPVAL,LAST,MPN ! 20 FORMAT(' PERMUTATION P-VALUE FROM TWO SAMPLE T-TEST' & & /' WITH SAMPLES M=',I3,3X,'N=',I3,6X,'P-VALUE',F8.4) ! ! LET'S USE THESE NUMBERS FOR A TEST PROBLEM DATA M/ 10/, N/ 5/ DATA X/ 52.43, 73.77, 74.88, 59.92, 47.15, 55.97, 59.54, 98.65, & & 52.27, 101.34, 10*0. / DATA Y / 74.69, 64.62, 49.47, 49.72, 40.58, 15*0. / ! OPEN( UNIT=6, FILE='t2perm.out' ) ! ! EVALUATE THE TEST STATISTIC FROM THE ORIGINAL DATA CALL TWOT(X,M,Y,N,TORIG) ! ! COMBINE THE TWO SAMPLES DO I = 1,M BOTH(I) = X(I) END DO ! LOOP ON I DO I = 1,N BOTH(I+M) = Y(I) END DO ! LOOP ON I ! MPN = M + N ! INITIALIZE COUNTER FOR P-VALUE KPVAL = 0 ! INITIALIZE SUBSET AND FLAG (LAST) DO I = 1,M LIST(I) = I END DO ! LOOP ON I LAST = 1 ! RETURN HERE TO START PROCESSING DO ! *** NOTE UNRESTRICTED DO *** EXIT AT END ! TWO INITIALIZATIONS TO MAKE SPLITTING EASIER LIST(M+1) = 0 KL = 1 ! BREAK COMBINED SAMPLE ACCORDING TO SUBSET IN LIST DO I = 1,MPN IF( I .EQ. LIST(KL) ) THEN X(KL) = BOTH(I) KL = KL + 1 ELSE Y(I+1-KL) = BOTH(I) END IF ! ( I .EQ. LIST(KL) ) END DO ! LOOP ON I ! EVALUATE T STATISTIC FOR THIS SUBSET CALL TWOT(X,M,Y,N,T) ! COUNT FOR P-VALUE IF( T .LT. TORIG ) KPVAL = KPVAL + 1 ! GET NEXT SUBSET CALL NXTKON(LIST,M,MPN,LAST) ! ARE WE DONE? IF( LAST .EQ. 1 ) EXIT END DO ! ! COMPUTE P-VALUE T = REAL(KPVAL) / RBICOF(MPN,N) ! WRITE OUT RESULTS WRITE(6,20) M,N,T ! STOP END PROGRAM T2PERM SUBROUTINE TWOT(X,M,Y,N,T) ! COMPUTE TWO-SAMPLE T-TEST ! ! X REAL VECTOR FIRST SAMPLE OF M OBSERVATIONS ! M INTEGER NUMBER OF OBSERVATIONS IN FIRST SAMPLE ! Y REAL VECTOR SECOND SAMPLE OF N OBSERVATIONS ! N INTEGER NUMBER OF OBSERVATIONS IN SECOND SAMPLE ! T REAL VALUE OF TWO SAMPLE T STATISTIC ! IMPLICIT NONE INTEGER, INTENT(IN) :: M,N REAL, DIMENSION(M), INTENT(IN) :: X REAL, DIMENSION(N), INTENT(IN) :: Y REAL, INTENT(OUT) :: T REAL XM,XS,YM,YS INTEGER I ! COMPUTE SAMPLE MEAN AND SUM OF SQUARES FOR X'S XM = 0. XS = 0. DO I = 1,M XM = XM + X(I) IF( I .GT. 1 ) XS = XS + (XM-I*X(I))*(XM-I*X(I)) / REAL(I*(I-1)) END DO ! LOOP ON I ! COMPUTE SAMPLE MEAN AND SUM OF SQUARES FOR X'S YM = 0. YS = 0. DO I = 1,N YM = YM + Y(I) IF( I .GT. 1 ) YS = YS + (YM-I*Y(I))*(YM-I*Y(I)) / REAL(I*(I-1)) END DO ! LOOP ON I ! MEANS FROM SUMS XM = XM / REAL(M) YM = YM / REAL(N) ! COMPUTE T STATISTIC T = (XM-YM)/SQRT( (1./REAL(M)+1./REAL(N)) * (XS+YS)/(N+M-2) ) RETURN END SUBROUTINE TWOT SUBROUTINE NXTKON(LIST,K,N,LAST) ! PRODUCES THE NEXT SUBSET OF SIZE K OUT OF A SET OF N ELEMENTS ! ENUMERATES THESE IN SEQUENCE LAST .EQ. 1 INDICATES LAST ONE ! ! LIST INTEGER VECTOR OF INDICES OF LENGTH K ! K INTEGER LENGTH OF INDEX VECTOR (SIZE OF SUBSET) ! N INTEGER SIZE OF SET ! LAST INTEGER 1 INDICATES LAST ONE ; 0 INDICATES MORE TO COME ! ! TO BE USED IN THE FOLLOWING FORMAT ! ! LAST = 1 ! ! START OF LOOP ! DO ! *** NOTE UNRESTRICTED DO *** ! CALL NXTKON(LIST,K,N,LAST) ! ! (CODE FOR USING SUBSET IN LIST) ! ! IF( LAST .EQ. 1 ) EXIT ! BRANCH OUT ON LAST ONE ! ! END OF LOOP ! END DO ! ! ! REMAINING CODE ! ! J F MONAHAN (FEB, 1987) DEPT OF STATISTICS, NCSU, RALEIGH, NC, USA ! RECODED MARCH 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: K,N INTEGER, INTENT(IN OUT) :: LAST INTEGER, DIMENSION(K), INTENT(IN OUT) :: LIST INTEGER I,KP,LKPMAX ! IF( LAST .EQ. 1) THEN ! INITIALIZE FOR FIRST ONE DO I = 1,K LIST(I) = I END DO ! LOOP ON I ! BEWARE OF CASE WITH ONLY ONE SUBSET IF( K .EQ. N ) RETURN LAST = 0 RETURN END IF ! ( LAST .EQ. 1 ) ! NOW DO THE USUAL CASE ! ! START WITH THE LAST INDEX KP = K LKPMAX = N+1 DO WHILE( KP .GT. 1 ) ! CAN WE INCREMENT THIS INDEX BY ONE? IF( LIST(KP) + 1 .LT. LKPMAX ) EXIT ! NO ROOM HERE, BACK UP TO NEXT ONE AND TRY LKPMAX = LIST(KP) KP = KP - 1 END DO ! WHILE( KP .GT. 1 ) ! CAN ALWAYS INCREMENT FIRST INDEX ON LAST GO ROUND ! MOVED INDEX KP UP ONE, NOW FIX THE REST LIST(KP) = LIST(KP) + 1 ! IS THIS THE LAST ONE? IF( LIST(1) .EQ. N - K + 1 ) LAST = 1 IF( KP .EQ. K ) RETURN DO WHILE ( KP .LT. K ) KP = KP + 1 LIST(KP) = LIST(KP-1) + 1 END DO ! WHILE ( KP .LT. K ) RETURN END SUBROUTINE NXTKON REAL FUNCTION RBICOF(I,J) ! COMPUTES BINOMIAL COEFFICIENTS, USUALLY BY TABLE ! IMPLICIT NONE INTEGER, INTENT(IN) :: I,J INTEGER K,L,N INTEGER, DIMENSION(21) :: PUSH REAL, DIMENSION(122) :: BICFTB REAL(KIND=8) SLGAMK ! ! TEST WHETHER IT IS IN TABLE - ! REMEMBER THE RESULT WILL GO HAYWIRE IN INTEGERS IF TOO BIG ! DATA PUSH/1,2,4,6,8,11,14,18,22,27,32,38,44,51,58,66,74,83,92, & & 102,112/ DATA BICFTB/1.,1.,1.,1.,2.,1.,3.,1.,4.,6.,1.,5.,10.,1.,6.,15., & & 20.,1.,7.,21.,35.,1.,8.,28.,56.,70.,1.,9.,36.,84.,126.,1.,10., & & 45.,120.,210.,252.,1.,11.,55.,165.,330.,462.,1.,12.,66.,220., & & 495.,792.,924.,1.,13.,78.,286.,715.,1287.,1716.,1.,14.,91.,364.& &,1001.,2002.,3003.,3432.,1.,15.,105.,455.,1365.,3003.,5005.,6435.& &,1.,16.,120.,560.,1820.,4368.,8008.,11440.,12870.,1.,17.,136., & & 680.,2380.,6188.,12376.,19448.,24310., 1.,18.,153.,816.,3060., & & 8568.,18564.,31824.,43758.,48620., 1.,19.,171.,969.,3876.,11628.& &,27132.,50388.,75582.,92378., 1.,20.,190.,1140.,4845.,15504., & & 38760.,77520.,125970.,167960.,184756. / ! ! TAKE CARE OF SPECIAL CASES RBICOF = 0. IF( ( J .LT. 0 ) .OR. ( J .GT. I ) ) RETURN RBICOF = 1. IF( ( J .EQ. 0 ) .OR. ( J .EQ. I ) ) RETURN ! GENERAL CASE ! L = MIN(J,I-J) IF( I .LE. 20 ) THEN ! LOOK IT UP IN THE TABLE K = PUSH(I+1) ! FIND THE STARTING POINT FOR I+1 RBICOF = BICFTB(K+L) ! THEN GET THE L-TH ELEMENT ELSE IF( L .LE. 20 ) THEN ! FOR BIG VALUES COMPUTE BY MULTIPLICATION RBICOF = 1. DO N = 1,L RBICOF = RBICOF*REAL(I-N+1) / REAL(N) END DO ! LOOP ON N ELSE ! FOR LONG MULTIPLICATIONS, USE STIRLING VIA SLGAMK RBICOF = REAL( SLGAMK(I) - SLGAMK(J) - SLGAMK(I-K) ) RBICOF = EXP( RBICOF ) END IF ! ( L .LE. 20 ) END IF ! ( I .LE. 20 ) RETURN END FUNCTION RBICOF REAL(KIND=8) FUNCTION SLGAMK(K) ! PRODUCES NATURAL LOGARITHM OF K! ! FOR SMALL VALUES OF K -- USES TABLE ! FOR LARGER VALUES -- USES STIRLING'S APPROXIMATION ! ! J F MONAHAN (SEPT, 1990) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: K REAL(KIND=8), DIMENSION(36) :: TAB REAL(KIND=8) F,DLF REAL(KIND=8), PARAMETER :: CHL2PI = .918938533205D0 ! log(2pi)/2 ! DATA TAB/ 0.D0, 0.6931471806D0, 1.7917594692D0, 3.1780538303D0& &, 4.7874917428D0, 6.5792512120D0, 8.5251613611D0,10.6046029027D0& &,12.8018274801D0,15.1044125731D0,17.5023078459D0,19.9872144957D0& &,22.5521638531D0,25.1912211827D0,27.8992713838D0,30.6718601061D0& &,33.5050734501D0,36.3954452080D0,39.3398841872D0,42.3356164608D0& &,45.3801388985D0,48.4711813518D0,51.6066755678D0,54.7847293981D0& &,58.0036052230D0,61.2617017610D0,64.5575386270D0,67.8897431372D0& &,71.2570389672D0,74.6582363488D0,78.0922235533D0,81.5579594561D0& &,85.0544670176D0,88.5808275422D0,92.1361756037D0,95.7196945421D0/ ! SLGAMK = 0.D0 IF( (K.EQ.0) .OR. (K.EQ.1) ) RETURN IF( K .LE. 36 ) THEN SLGAMK = TAB(K) RETURN ! USE STIRLING'S FORMULA ELSE F = DBLE(K) DLF = LOG(F) SLGAMK = ( (F+.5D0)*DLF - F ) & & + ( CHL2PI + ( 1.D0 - 1.D0/(30.D0*F*F) )/(12.D0*F) ) END IF ! ( K .LE. 36 ) ! RETURN END FUNCTION SLGAMK ! *** end of filename t2perm.f95 *********************