! Last change: DOS 3 Aug 2000 6:10 pm ! *** copyright 2000 *** ! *** filename nxtkon.f95 *** John F. Monahan ** ! ********************** PROGRAM PNXTKON ! TEST PROGRAM FOR NXTKON ! ENUMERATES ALL SUBSETS OF SIZE K OUT OF N ELEMENTS ! ! COMPUTE THE NULL DISTRIBUTION OF THE WILCOXON TWO SAMPLE STATISTIC ! WITH SAMPLE SIZES N1 = K AND N2 = N-K ! IMPLICIT NONE INTEGER, DIMENSION(10) :: LIST INTEGER, DIMENSION(200) :: LISTPR INTEGER J,N,K,LAST,INDEX,MAX,S INTEGER KKOMBO ! 20 FORMAT(/' TEST OF NXTKON -- COMPUTES ALL (N CHOOSE K) SUBSETS' & & /' FOR K=',I3,' N=',I3,' AND COMPUTES NULL DISTRIBUTION'& & /' OF WILCOXON TWO SAMPLE STATISTIC WITH N1=K, N1+N2=N'/) 21 FORMAT(10I5) 22 FORMAT(' INDEX=',I5,3X,'LIST',10I4) ! OPEN( UNIT=6, FILE='nxtkon.out' ) ! USE THESE VALUES TO TEST N = 12 K = 4 ! INITIALIZE LIST OF PROBS (REALLY JUST COUNTS) MAX = ( N * (N+1) )/2 DO J = 1,MAX LISTPR(J) = 0 END DO ! LOOP ON J ! INITIALIZE SUBSET COUNTER LAST = 1 ! START OF LOOP DO ! *** NOTE UNRESTRICTED DO *** CALL NXTKON(LIST,K,N,LAST) ! ! COMPUTE VALUE OF TEST STATISTIC S = 0 DO J = 1,K S = S + LIST(J) END DO ! LOOP ON J ! INCREMENT THIS VALUE BY ONE LISTPR(S) = LISTPR(S) + 1 ! COMBINATION INDEX INDEX = KKOMBO(LIST,K,N) ! DEBUG WRITE WRITE(6,22) INDEX,(LIST(J),J=1,K) ! IF( LAST .EQ. 1 ) EXIT ! BRANCH OUT ON LAST ONE ! ! END OF LOOP END DO ! ! WRITE OUT DISTRIBUTION WRITE(6,20) K,N WRITE(6,21) (LISTPR(J),J=1,MAX) STOP END PROGRAM PNXTKON 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 INTEGER FUNCTION KKOMBO(LIST,K,N) ! PRODUCES A UNIQUE INTEGER FOR EACH SUBSET OF SIZE K OUT OF N ! DESCRIBED IN THE VECTOR LIST *** MODIFICATION OF KNOTT'S ! KKOMBO RANGES IN VALUE FROM 1 FOR ( 1, 2, ..., K ) ! TO ! (N CHOOSE K) FOR ( N-K+1, ..., N-1, N ) ! ! LIST INTEGER VECTOR OF INDICES OF LENGTH K ! K INTEGER LENGTH OF INDEX VECTOR (SIZE OF SUBSET) ! N INTEGER SIZE OF SET ! ! SEE G D KNOTT (1974) "A NUMBERING SYSTEM FOR COMBINATIONS," ! CACM, VOLUME 17, PP. 45-46 ! ! J F MONAHAN (AUGUST 1998) DEPT OF STATISTICS, NC STATE UNIV ! RECODED MARCH 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: K,N INTEGER, DIMENSION(K), INTENT(IN) :: LIST INTEGER I,LISTM1 INTEGER IBICOF ! INITIALIZE KKOMBO = IBICOF(N,K) ! COMPUTE KNOTT'S FUNCTION FOR REVERSE COMBO DO I = 1,K LISTM1 = N - LIST(K+1-I) ! AND SUBTRACT FROM MAX KKOMBO = KKOMBO - IBICOF(LISTM1,I) END DO ! LOOP ON I RETURN END FUNCTION KKOMBO INTEGER FUNCTION IBICOF(I,J) ! COMPUTES BINOMIAL COEFFICIENTS, USUALLY BY TABLE ! INTEGER, INTENT(IN) :: I,J INTEGER K,L,N INTEGER, DIMENSION(21) :: PUSH ! TABLE OF POINTERS INTEGER, DIMENSION(122) :: BICFTB ! ! 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 IBICOF = 0 IF( ( J .LT. 0 ) .OR. ( J .GT. I ) ) RETURN IBICOF = 1 IF( ( J .EQ. 0 ) .OR. ( J .EQ. I ) ) RETURN ! GENERAL CASE L = MIN(J,I-J) IF( I .LE. 20 ) THEN ! TABLE LOOKUP K = PUSH(I+1) ! FIND START POINT FOR I IBICOF = BICFTB(K+L) ! PULL NUMBER FROM TABLE ELSE ! MULTIPLY ! FOR BIG VALUES COMPUTE BY MULTIPLICATION IBICOF = 1 DO N = 1,L IBICOF = (IBICOF*(I-N+1)) / N END DO ! LOOP ON N END IF ! ( I .LE. 20 ) RETURN END FUNCTION IBICOF ! *** end of filename nxtkon.f95 *********************