PARAMETER (N=30,M=20) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Z(N,M),X(N,M),R(M,M),RR(M,M),AV(M),SD(M),XX(N),U(N),V(N) OPEN(7,FILE='INDICES.TXT') DO I=1,N READ(7,*)(Z(I,J),J=1,M) WRITE(*,*)I,(Z(I,J),J=1,M) ENDDO CLOSE(7) OPEN(7,FILE='CORIND.TXT') WRITE(*,*)'FEED CHOICE OF CORRELATION. ZERO (0) IS KARL PEARSON R' WRITE(*,*)'(1): RANK CORRELATION, (2): SIGNUM CORRELATION' WRITE(*,*)'(3): BRADLEY CORRELATION, (4) SHEVLYAKOV CORRELATION' WRITE(*,*)'(5): CAMPBELL CORRELATION MATRIX' WRITE(*,*)'(6): CAMPBELL CORRELATION MATRIX (MEDIAN-BASED WEIGHT)' READ(*,*) NTYPE DO J=1,M DO I=1,N X(I,J)=Z(I,J) ENDDO ENDDO C ----------------------------------------------------------------- IF(NTYPE.EQ.0) THEN WRITE(*,*) 'KARL PEARSON CORRELATION MATRIX' CALL PROD(X,N,M,AV,SD,R) ENDIF C ----------------------------------------------------------------- IF(NTYPE.EQ.1) THEN WRITE(*,*) 'SPEARMAN RANK CORRELATION MATRIX' DO J=1,M DO I=1,N XX(I)=Z(I,J) ENDDO CALL RANK(XX,N) ! RANK TRANSFORMATION OF X DO I=1,N X(I,J)=XX(I) ENDDO ENDDO CALL PROD(X,N,M,AV,SD,R) ENDIF C ------------------------------------------------------------------ IF(NTYPE.EQ.2) THEN WRITE(*,*) 'SIGNUM CORRELATION MATRIX' DO J=1,M DO I=1,N XX(I)=Z(I,J) ENDDO CALL MEDIAN(XX,N,AMED,AMDEV) DO I=1,N X(I,J)=1.D0 DEV=(Z(I,J)-AMED) IF(DABS(DEV).GT.1.D-06) X(I,J)=DEV/DABS(DEV) ENDDO ENDDO DO J=1,M DO JJ=1,M R(J,JJ)=0.D0 DO I=1,N CONC=X(I,J)*X(I,JJ) C IF(CONC.EQ.1) R(J,JJ)=R(J,JJ)+CONC R(J,JJ)=R(J,JJ)+CONC ENDDO R(J,JJ)=R(J,JJ)/N ENDDO ENDDO DO J=1,M WRITE(7,*)(R(J,JJ),JJ=1,M) ENDDO WRITE(*,*)'------- SIGNUM CORRELATION MATRIX (PROPER) -----------' CALL PROD(X,N,M,AV,SD,R) WRITE(*,*)'OVER' ENDIF C ----------------------------------------------------------------- IF(NTYPE.EQ.3) THEN WRITE(*,*) 'BRADLEY CORRELATION MATRIX' DO J=1,M DO I=1,N XX(I)=Z(I,J) ENDDO CALL MEDIAN(XX,N,AMED,AMDEV) DO I=1,N X(I,J)=0.D0 DEV=(Z(I,J)-AMED) X(I,J)=DEV/AMDEV ENDDO ENDDO DO J=1,M DO JJ=1,M R(J,JJ)=0.D0 S1=0.D0 S2=0.D0 DO I=1,1,N S1=S1+DABS(X(I,J)+X(I,JJ)) - DABS(X(I,J)-X(I,JJ)) S2=S2+DABS(X(I,J))+ DABS(X(I,JJ)) ENDDO R(J,JJ)=S1/S2 ENDDO ENDDO DO J=1,M WRITE(7,*)(R(J,JJ),JJ=1,M) ENDDO WRITE(*,*)'------------------------------------------------------' ENDIF C ----------------------------------------------------------------- IF(NTYPE.EQ.4) THEN WRITE(*,*) 'SHEVLYAKOV CORRELATION MATRIX' DO J=1,M DO I=1,N XX(I)=Z(I,J) ENDDO CALL MEDIAN(XX,N,AMED,AMDEV) DO I=1,N XX(I)=DABS(Z(I,J)-AMED) ENDDO CALL MEDIAN(XX,N,HMED,HMDEV) DO I=1,N X(I,J)=(Z(I,J)-AMED)/HMED ENDDO ENDDO DO J=1,M DO JJ=1,M R(J,JJ)=0.D0 DO I=1,N U(I)=DABS(X(I,J)+X(I,JJ)) V(I)=DABS(X(I,J)-X(I,JJ)) ENDDO CALL MEDIAN(U,N,UMED,UMDEV) CALL MEDIAN(V,N,VMED,VMDEV) R(J,JJ)=(UMED**2-VMED**2)/(UMED**2+VMED**2) ENDDO ENDDO DO J=1,M WRITE(7,*)(R(J,JJ),JJ=1,M) ENDDO WRITE(*,*)'------------------------------------------------------' ENDIF C ----------------------------------------------------------------- IF(NTYPE.EQ.5) THEN WRITE(*,*) 'CAMPBELL CORRELATION MATRIX' DO I=1,N DO J=1,M X(I,J)=Z(I,J) ENDDO ENDDO CALL RCAMPBELL(X,N,M,R) ENDIF C ----------------------------------------------------------------- IF(NTYPE.EQ.6) THEN WRITE(*,*) 'CAMPBELL CORRELATION MATRIX (MEDIAN-BASED)' DO I=1,N DO J=1,M X(I,J)=Z(I,J) ENDDO ENDDO CALL MCAMPBELL(X,N,M,R) ENDIF C ----------------------------------------------------------------- 1 FORMAT(10F7.3) CLOSE(7) END C ----------------------------------------------------------------- SUBROUTINE PROD(X,N,M,AV,SD,R) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N,M),R(M,M),AV(M),SD(M) DO J=1,M AV(J)=0.D0 SD(J)=0.D0 DO I=1,N AV(J)=AV(J)+X(I,J) SD(J)=SD(J)+X(I,J)**2 ENDDO SD(J)=DSQRT((N*SD(J)-AV(J)*AV(J))/N**2) AV(J)=AV(J)/N ENDDO DO J=1,M DO JJ=1,M R(J,JJ)=0.D0 DO I=1,N R(J,JJ)=R(J,JJ)+X(I,J)*X(I,JJ) ENDDO R(J,JJ)=(R(J,JJ)/N-AV(J)*AV(JJ))/(SD(J)*SD(JJ)) ENDDO ENDDO DO J=1,M WRITE(7,*)(R(J,JJ),JJ=1,M) ENDDO WRITE(*,1)(AV(J),J=1,M) WRITE(*,1)(SD(J),J=1,M) 1 FORMAT(8F9.5) WRITE(*,*)'------------------------------------------------------' RETURN END C ----------------------------------------------------------------- SUBROUTINE RANK(X,N) PARAMETER (NMAX=1000) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N), SL(NMAX) DO I=1,N SL(I)=DFLOAT(I) ENDDO DO I=1,N-1 DO II=I+1,N IF(X(I).LT.X(II)) THEN T=X(I) X(I)=X(II) X(II)=T T=SL(I) SL(I)=SL(II) SL(II)=T ENDIF ENDDO ENDDO DO I=1,N X(I)=I ENDDO DO I=1,N-1 DO II=I+1,N IF(SL(I).GT.SL(II)) THEN T=X(I) X(I)=X(II) X(II)=T T=SL(I) SL(I)=SL(II) SL(II)=T ENDIF ENDDO ENDDO RETURN END C ----------------------------------------------------------------- SUBROUTINE MEDIAN(X,N,A,V) ! ------------------------------------ C SUBROUTINE MEDIAN : FINDS MEDIAN (A) AND MEAN DEVIATION (V) OF A C GIVEN VARIATE, VARIATE X(N) PARAMETER (NMAX=1000) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),Z(NMAX) C STORE X IN Z DO I=1,N Z(I)=X(I) ENDDO C ARRANGE Z IN AN ASCENDING ORDER DO I=1,N-1 DO J=I+1,N IF(Z(I).GT.Z(J)) THEN ! EXCHANGE TEMP=Z(I) Z(I)=Z(J) Z(J)=TEMP ENDIF ENDDO ENDDO K=(N+1)/2 ! K IS OBTAINED AS INT((N+1)/2.0D0) A=(Z(K)+Z(N+1-K))/2.D0 ! GIVES MEDIAN FOR ODD AS WELL AS EVEN N C FIND MEAN DEVIATION V=0.D0 DO I=1,N V=V+DABS(Z(I)-A) ! A IS MEDIAN ENDDO V=V/N ! V IS MEAN DEVIATION FROM MEDIAN C WRITE(*,*)'MEDIAN =',A,' MEAN DEVIATION =',V RETURN END C ------------------------------------------------------------------ C CAMPBELL CORRELATION MATRIX SUBROUTINE RCAMPBELL(X,N,M,V) PARAMETER(NV=30,MV=20,ITRN=200) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N,M),V(M,M),AV(MV),W(NV),XD(MV) DIMENSION D(NV),VV(MV,MV),DN(NV) DATA B1,B2/2, 1.5/ C SOME DEFINITIONS D0=DSQRT(DFLOAT(M))+B1/DSQRT(2.D0) B22=B2**2 NSKIP=1 IF(NSKIP.NE.1) THEN ! DO NOT STANDARDIZE THE VARIABLES C STANDARDIZE DO J=1,M AV(J)=0.D0 XD(J)=0.D0 DO I=1,N AV(J)=AV(J)+X(I,J) XD(J)=XD(J)+X(I,J)**2 ENDDO AV(J)=AV(J)/N XD(J)=DSQRT(XD(J)/N-AV(J)**2) ENDDO DO J=1,M DO I=1,N X(I,J)=(X(I,J)-AV(J))/XD(J) ENDDO ENDDO ENDIF C INITIALIZE WEIGHT VECTOR BY UNITY DO I=1,N W(I)=1.D0 ENDDO C FIND SUM OF WEIGHTS DO ITER=1,ITRN SW=0.D0 SSW=0.D0 DO I=1,N SW=SW+W(I) SSW=SSW+W(I)**2 ENDDO SSW=SSW-1.D0 C COMPUTE MEAN VECTOR AND COVARIANCE MATRIX DO J=1,M AV(J)=0.D0 DO I=1,N AV(J)=AV(J)+X(I,J)*W(I) ENDDO AV(J)=AV(J)/SW ENDDO DO J=1,M DO JJ=J,M V(J,JJ)=0.D0 DO I=1,N V(J,JJ)=V(J,JJ)+(X(I,J)-AV(J))*(X(I,JJ)-AV(JJ))*W(I)**2 ENDDO V(J,JJ)=V(J,JJ)/SSW IF(J.NE.JJ) V(JJ,J)=V(J,JJ) ENDDO ENDDO DO J=1,M DO JJ=1,M VV(J,JJ)=V(J,JJ) ENDDO ENDDO C INVERT V CALL MINV(V,M,DD) ! ON RETURN V IS INVERTED V DO I=1,N D(I)=0.D0 DO J=1,M XD(J)=0.D0 DO JJ=1,M XD(J)=XD(J)+(X(I,JJ)-AV(JJ))*V(JJ,J) ENDDO ENDDO DD=0.D0 DO J=1,M DD=DD+XD(J)*(X(I,J)-AV(J)) ENDDO DD=DSQRT(DD) D(I)=DD DN(I)=DD ENDDO DO I=1,N IF(D(I).LE.D0)THEN WD= D(I) ELSE WD=D0*DEXP(-0.5D0*(D(I)-D0)**2/B22) ENDIF W(I)=1.D0 IF(DABS(D(I)).GT.1.0D-05) W(I)=WD/D(I) ENDDO ENDDO DO J=1,M DO JJ=1,M V(J,JJ)=VV(J,JJ)/DSQRT(VV(J,J)*VV(JJ,JJ)) ENDDO ENDDO DO J=1,M WRITE(7,*)(V(J,JJ),JJ=1,M) ENDDO 1 FORMAT(8F9.5) WRITE(*,*)'------------------------------------------------------' RETURN END C ----------------------------------------------------------------- C ------------------------------------------------------------------ C CAMPBELL CORRELATION MATRIX SUBROUTINE MCAMPBELL(X,N,M,V) PARAMETER(NV=30,MV=20,ITRN=200) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N,M),V(M,M),AV(MV),W(NV),XD(MV) DIMENSION D(NV),VV(MV,MV),DN(NV) DATA B1,B2/2, 1.5/ C SOME DEFINITIONS D0=DSQRT(DFLOAT(M))+B1/DSQRT(2.D0) B22=B2**2 NSKIP=1 IF(NSKIP.NE.1) THEN ! DO NOT STANDARDIZE THE VARIABLES C STANDARDIZE DO J=1,M AV(J)=0.D0 XD(J)=0.D0 DO I=1,N AV(J)=AV(J)+X(I,J) XD(J)=XD(J)+X(I,J)**2 ENDDO AV(J)=AV(J)/N XD(J)=DSQRT(XD(J)/N-AV(J)**2) ENDDO DO J=1,M DO I=1,N X(I,J)=(X(I,J)-AV(J))/XD(J) ENDDO ENDDO ENDIF C INITIALIZE WEIGHT VECTOR BY UNITY DO I=1,N W(I)=1.D0 ENDDO C FIND SUM OF WEIGHTS DO ITER=1,ITRN SW=0.D0 SSW=0.D0 DO I=1,N SW=SW+W(I) SSW=SSW+W(I)**2 ENDDO SSW=SSW-1.D0 C COMPUTE MEAN VECTOR AND COVARIANCE MATRIX DO J=1,M AV(J)=0.D0 DO I=1,N AV(J)=AV(J)+X(I,J)*W(I) ENDDO AV(J)=AV(J)/SW ENDDO DO J=1,M DO JJ=J,M V(J,JJ)=0.D0 DO I=1,N V(J,JJ)=V(J,JJ)+(X(I,J)-AV(J))*(X(I,JJ)-AV(JJ))*W(I)**2 ENDDO V(J,JJ)=V(J,JJ)/SSW IF(J.NE.JJ) V(JJ,J)=V(J,JJ) ENDDO ENDDO DO J=1,M DO JJ=1,M VV(J,JJ)=V(J,JJ) ENDDO ENDDO C INVERT V CALL MINV(V,M,DD) ! ON RETURN V IS INVERTED V DO I=1,N D(I)=0.D0 DO J=1,M XD(J)=0.D0 DO JJ=1,M XD(J)=XD(J)+(X(I,JJ)-AV(JJ))*V(JJ,J) ENDDO ENDDO DD=0.D0 DO J=1,M DD=DD+XD(J)*(X(I,J)-AV(J)) ENDDO DD=DSQRT(DD) D(I)=DD DN(I)=DD ENDDO CALL MEDIAN(DN,N,DNA,DNV) DO I=1,N DN(I)=DABS(DN(I)-DNA) ENDDO CALL MEDIAN(DN,N,DNAA,DNVV) DNAA=DNAA/0.6745 DO I=1,N W(I)=0.D0 DX=DABS(D(I)-DNA) IF(DX.LE.DNAA) W(I)=1.D0 IF(DX.LE.2*DNAA.AND.DX.GT.DNAA) W(I)=.25D0 IF(DX.LE.3*DNAA.AND.DX.GT.2*DNAA) W(I)=0.1111111D0 c IF(DX.LE.4*DNAA.AND.DX.GT.3*DNAA) W(I)=0.0625D0 ENDDO ENDDO DO J=1,M DO JJ=1,M V(J,JJ)=VV(J,JJ)/DSQRT(VV(J,J)*VV(JJ,JJ)) ENDDO ENDDO DO J=1,M WRITE(7,*)(V(J,JJ),JJ=1,M) ENDDO 1 FORMAT(8F9.5) WRITE(*,*)'------------------------------------------------------' RETURN END C ----------------------------------------------------------------- C SUBROUTINE FOR MATRIX INVERSION SUBROUTINE MINV(A,N,D) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(N,N) U=1.D0 D=U DO I=1,N D=D*A(I,I) A(I,I)=U/A(I,I) DO J=1,N IF(I.NE.J) A(J,I)=A(J,I)*A(I,I) ENDDO DO J=1,N DO K=1,N IF(I.NE.J.AND.K.NE.I) A(J,K)=A(J,K)-A(J,I)*A(I,K) ENDDO ENDDO DO J=1,N IF(J.NE.I) A(I,J)= -A(I,J)*A(I,I) ENDDO ENDDO RETURN END C -----------------------------------------------------------------