C @(#)rankcorr.for 17.1.1.1 (ESO-DMD) 01/25/02 17:12:48 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 1 JUL 1989 C C.LANGUAGE: F77+ESOext C C.AUTHOR: M.PERON C C.IDENTIFICATION C C Program RANKCORR C C.KEYWORDS C C KENDALL, SPEARMANN COEFFICIENT C RANK CORRELATION C.PURPOSE C calculates the Spearman and Kendall rank correlation coefficients and C the significance level of the test C C C 010621 last modif C C----------------------------------------------------------- C PROGRAM RANKCORR C IMPLICIT NONE C REAL D,PROB,PROBD,PROBRS,TAU,Z,ZD,RS C INTEGER*8 ADDRIN(4) C INTEGER MADRID INTEGER PARVAL,NPAR,ISTAT,TID1,NCOL,NROW1,NACOL,NAROW,COL1,COL2 INTEGER IAV,NREAL INTEGER NSC,UNI,NULO,NBYTES,KUN C CHARACTER*80 INTAB CHARACTER*17 COLIND,COLDEP CHARACTER*80 OUTPUT1,OUTPUT2,OUTPUT3,OUT,OUTPUT0,OUTPUT4 CHARACTER*1 ACTION C INCLUDE 'MID_INCLUDE:TABLES.INC' COMMON/VMR/MADRID(1) C INCLUDE 'MID_INCLUDE:TABLED.INC' DATA PARVAL/8/, TID1 /-1/ DATA OUTPUT0 1 /' Null Hypothesis: The two data sets are not associated'/ DATA OUTPUT1 1 /' Kendall Rank-Order coeff '/ DATA OUTPUT2 1 /' Spearman Rank-Order coeff '/ DATA OUTPUT3 2 /' ---------------------------------------------------------'/ DATA OUTPUT4 3 /'Probability of exceeding this value under the 4 Null Hypothesis'/ C CALL STSPRO('RANKCORR') c get the parameters CALL TDPGET(PARVAL,NPAR,ISTAT) IF (ISTAT.NE.0) GOTO 80 INTAB = TPARBF(1) COLIND = TPARBF(2) COLDEP = TPARBF(3) c open the table CALL TBTOPN(INTAB,F_I_MODE,TID1,ISTAT) IF(ISTAT.NE.0) GOTO 80 CALL TBIGET(TID1,NCOL,NROW1,NSC,NACOL,NAROW,ISTAT) CALL TBCSER(TID1,COLIND,COL1,ISTAT) IF (COL1.LT.0) THEN CALL STTPUT('Column not found...',ISTAT) GOTO 80 ENDIF CALL TBCSER(TID1,COLDEP,COL2,ISTAT) IF (COL2.LT.0) THEN CALL STTPUT('Column not found...',ISTAT) GOTO 80 ENDIF c Choose between Kendall and Spearman CALL STKRDC('ACTION',1,1,1,IAV,ACTION,UNI,NULO,ISTAT) CALL UPCAS(ACTION,ACTION) IF (ACTION.EQ.'K')THEN !KEndall coeff NBYTES = 4*NROW1 CALL TDMGET(NBYTES,ADDRIN(1),ISTAT) CALL TDMGET(NBYTES,ADDRIN(2),ISTAT) CALL TMAP2(TID1,NROW1,COL1,COL2,MADRID(ADDRIN(1)), 1 MADRID(ADDRIN(2)),NREAL) c calculates the coefficient and the proba CALL KEND(MADRID(ADDRIN(1)),MADRID(ADDRIN(2)),NREAL,TAU,Z,PROB) WRITE(OUT,100)TAU CALL STTPUT(OUTPUT0,ISTAT) CALL STTPUT(OUTPUT3,ISTAT) CALL STTPUT(OUTPUT1,ISTAT) CALL STTPUT(OUT,ISTAT) CALL STTPUT(OUTPUT3,ISTAT) CALL STTPUT(OUTPUT4,ISTAT) CALL STTPUT(OUTPUT3,ISTAT) WRITE(OUT,100)PROB CALL STTPUT(OUT,ISTAT) c write results in keyword CALL STKWRR('OUTPUTR',TAU,1,1,KUN,ISTAT) CALL STKWRR('OUTPUTR',Z,2,1,KUN,ISTAT) CALL STKWRR('OUTPUTR',PROB,3,1,KUN,ISTAT) ELSE !Spearmann coeff NBYTES = 4*NROW1 CALL TDMGET(NBYTES,ADDRIN(1),ISTAT) CALL TDMGET(NBYTES,ADDRIN(2),ISTAT) CALL TDMGET(NBYTES,ADDRIN(3),ISTAT) CALL TDMGET(NBYTES,ADDRIN(4),ISTAT) CALL TMAP2(TID1,NROW1,COL1,COL2,MADRID(ADDRIN(1)), 1 MADRID(ADDRIN(2)),NREAL) c calculates the coefficient and the proba CALL SPEAR(MADRID(ADDRIN(1)),MADRID(ADDRIN(2)),NREAL, 1 MADRID(ADDRIN(3)),MADRID(ADDRIN(4)),D,ZD,PROBD,RS,PROBRS) WRITE(OUT,100)RS CALL STTPUT(OUTPUT0,ISTAT) CALL STTPUT(OUTPUT3,ISTAT) CALL STTPUT(OUTPUT2,ISTAT) CALL STTPUT(OUT,ISTAT) CALL STTPUT(OUTPUT3,ISTAT) CALL STTPUT(OUTPUT4,ISTAT) WRITE(OUT,100)PROBRS CALL STTPUT(OUT,ISTAT) CALL STTPUT(OUTPUT3,ISTAT) c write results in keyword CALL STKWRR('OUTPUTR',RS,1,1,KUN,ISTAT) CALL STKWRR('OUTPUTR',ZD,2,1,KUN,ISTAT) CALL STKWRR('OUTPUTR',PROBRS,3,1,KUN,ISTAT) ENDIF CALL TBTCLO(TID1,ISTAT) C That's all CALL STSEPI 80 CONTINUE 100 FORMAT(1X,G15.7) END SUBROUTINE KEND(X,Y,N,TAU,Z,PROB) INTEGER IS,J,K,N1,N2,N REAL X(N),Y(N),TAU,Z,PROB,A1,A2,AA REAL VAR,ERFCC N1 = 0 N2 = 0 IS = 0 DO 10 J = 1,N-1 DO 20 K = J+1,N A1 = X(J)-X(K) A2 = Y(J)-Y(K) AA = A1*A2 IF (AA.NE.0)THEN N1 = N1+1 N2 = N2+1 IF (AA.GT.0.)THEN IS = IS+1 ELSE IS = IS-1 ENDIF ELSE IF (A1.NE.0.) N1 = N1+1 IF (A2.NE.0.) N2 = N2+1 ENDIF 20 CONTINUE 10 CONTINUE TAU = FLOAT(IS)/SQRT(FLOAT(N1)*FLOAT(N2)) VAR = (4.*N+10.)/(9.*N*(N-1.)) Z = TAU/SQRT(VAR) PROB = ERFCC(ABS(Z/SQRT(2.))) END SUBROUTINE SPEAR(X,Y,N,WX,WY,D,ZD,PROBD,RS,PROBRS) INTEGER J,N REAL X(N),Y(N),WX(N),WY(N),D,ZD,PROBD,PROBRS,SF,SG,IRANKX(100) REAL IRANKY(100),RS,AVED,EN,EN3N,X2,Y2,VARD,BETAI,ERF REAL DF,T,FAC,RS1 DO 10 J=1,N WX(J) = X(J) WY(J) = Y(J) 10 CONTINUE CALL INDEXX(N,X,WX) CALL CRANK(N,X,WX,IRANKX,SF) CALL INDEXX(N,Y,WY) CALL CRANK(N,Y,WY,IRANKY,SG) D=0. DO 20 J=1,N D = D+(IRANKX(J)-IRANKY(J))**2 20 CONTINUE EN = N EN3N = EN**3-EN AVED = EN3N/6.-(SF+SG)/12. FAC = (1.-SF/EN3N)*(1.-SG/EN3N) X2 = EN3N/12.-SF/12. Y2 = EN3N/12 -SG/12. RS = 0.5*(AVED-D)/SQRT(X2*Y2) VARD = ((EN-1.)*EN**2*(EN+1.)**2/36.)*FAC ZD = (D-AVED)/SQRT(VARD) PROBD = ERF(ABS(ZD)) RS1 = (1.-(6./EN3N)*(D+0.5*(SF+SG)))/FAC T = RS*SQRT((EN-2.)/((1.+RS)*(1.-RS))) DF = EN-2 PROBRS = BETAI(0.5*DF,0.5,DF/(DF+T**2)) RETURN END