C @(#)prc.for 17.1.1.1 (ES0-DMD) 01/25/02 17:55:51 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 CC CC *STATIS* : STANDARDIZATION OF THE OBSERVATIONS AND CC CALCULATION OF SOME STATISTICS CC SUBROUTINE STATIS(X,Y,XMED,XMAD,RESDU,WEIGHTS,AW, 1 NMVAL,C,H,JNDVC,NSTOP,NVAR,NVAD,NVSB,NCAS,JCST,JPRT, 1 NDXX,NDXY,LUB,NMXV,PREC,LAB) INCLUDE 'MID_REL_INCL:implicit.inc' INTEGER J,JCST,JNC,JPRT,JS,JSM,JSPA,JSPB,JSPP,LUB INTEGER NCAS,NDXX,NDXY,NMXV,NVAD,NSTOP,NVAR,NVSB REAL AL,AMDAN,DIAG,SPEAR,PEARC,PREC C REAL X(NDXX,NDXY),Y(NDXY),XMED(NDXX),XMAD(NDXX) REAL RESDU(NDXY),WEIGHTS(NDXY),AW(NDXY),C(NMXV,NDXX) INTEGER NMVAL(NDXY),JNDVC(NDXX) C DIMENSION X(NDXX,NDXY),Y(NDXY),XMED(NDXX),XMAD(NDXX) C DIMENSION RESDU(NDXY),WEIGHTS(NDXY),AW(NDXY),NMVAL(NDXY) C DIMENSION C(NMXV,NDXX),JNDVC(NDXX) DOUBLE PRECISION H(NMXV,NDXX) CHARACTER*10 LAB(NDXX) CC-----IN THIS SUBROUTINE THE VECTOR "RESDU" WILL BE USED AS A CC-----WORKING VECTOR AL=NCAS IF (JCST.NE.0) GOTO 60 DO 50 J=1,NVAD XMED(J)=0.0 DO 10 JNC=1,NCAS 10 RESDU(JNC)=ABS(X(J,JNC)) XMAD(J)=AMDAN(AW,NDXY,RESDU,NCAS)*1.4826 IF (ABS(XMAD(J)).GT.PREC) GOTO 30 XMAD(J)=0.0 DO 20 JNC=1,NCAS 20 XMAD(J)=XMAD(J)+RESDU(JNC) XMAD(J)=(XMAD(J)/AL)*1.2533 IF (ABS(XMAD(J)).GT.PREC) GOTO 30 WRITE(LUB,8000) LAB(J) NSTOP=1 RETURN 30 DO 40 JNC=1,NCAS X(J,JNC)=X(J,JNC)/XMAD(J) 40 CONTINUE 50 CONTINUE WRITE(LUB,8015) GOTO 150 60 DO 120 J=1,NVAD IF(J.EQ.NVAR) GOTO 120 DO 70 JNC=1,NCAS 70 RESDU(JNC)=X(J,JNC) XMED(J)=AMDAN(AW,NDXY,RESDU,NCAS) DO 80 JNC=1,NCAS 80 RESDU(JNC)=ABS(RESDU(JNC)-XMED(J)) XMAD(J)=AMDAN(AW,NDXY,RESDU,NCAS)*1.4826 IF (ABS(XMAD(J)).GT.PREC) GOTO 100 XMAD(J)=0.0 DO 90 JNC=1,NCAS 90 XMAD(J)=XMAD(J)+RESDU(JNC) XMAD(J)=(XMAD(J)/AL)*1.2533 IF (ABS(XMAD(J)).GT.PREC) GOTO 100 WRITE(LUB,8000) LAB(J) NSTOP=1 RETURN 100 DO 110 JNC=1,NCAS X(J,JNC)=(X(J,JNC)-XMED(J))/XMAD(J) 110 CONTINUE 120 CONTINUE XMED(NVAR)=1.0 WRITE(LUB,8010) IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NVSB),LAB(NVAD) IF (JCST.EQ.0) WRITE(LUB,8020) (LAB(J),J=1,NVAD) IF ((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1)) 1 GOTO 130 IF (JCST.EQ.0) WRITE(LUB,8030)(XMED(J),J=1,NVAD) IF (JCST.EQ.1) WRITE(LUB,8030)(XMED(J),J=1,NVSB),XMED(NVAD) GOTO 140 130 WRITE(LUB,8030)(XMED(J),J=1,6) IF (JCST.EQ.0) WRITE(LUB,8040)(XMED(J),J=7,NVAD) IF (JCST.EQ.1.AND.NVSB.GE.7) 1 WRITE(LUB,8040)(XMED(J),J=7,NVSB),XMED(NVAD) IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8040) XMED(NVAD) 140 XMED(NVAR)=0.0 WRITE(LUB,8050) IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NVSB),LAB(NVAD) IF (JCST.EQ.0) WRITE(LUB,8020) (LAB(J),J=1,NVAD) 150 IF((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1)) 1 GOTO 160 IF (JCST.EQ.0) WRITE(LUB,8030)(XMAD(J),J=1,NVAD) IF (JCST.EQ.1) WRITE(LUB,8030)(XMAD(J),J=1,NVSB),XMAD(NVAD) GOTO 170 160 WRITE(LUB,8030)(XMAD(J),J=1,6) IF (JCST.EQ.0) WRITE(LUB,8040)(XMAD(J),J=7,NVAD) IF (JCST.EQ.1.AND.NVSB.GE.7) 1 WRITE(LUB,8040)(XMAD(J),J=7,NVSB),XMAD(NVAD) IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8040) XMAD(NVAD) 170 IF(JCST.NE.0) XMAD(NVAR)=1.0 IF (JPRT.LT.2) GOTO 200 WRITE(LUB,8060) IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NVSB),LAB(NVAD) IF (JCST.EQ.0) WRITE(LUB,8020) (LAB(J),J=1,NVAD) DO 190 JNC=1,NCAS IF((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1)) 1 GOTO 180 IF (JCST.EQ.1) WRITE(LUB,8070) NMVAL(JNC),(X(J,JNC),J=1,NVSB), 1 X(NVAD,JNC) IF (JCST.EQ.0) WRITE(LUB,8070) NMVAL(JNC),(X(J,JNC),J=1,NVAD) GOTO 190 180 WRITE(LUB,8070) NMVAL(JNC),(X(J,JNC),J=1,6) IF (JCST.EQ.1.AND.NVSB.GE.7) 1 WRITE(LUB,8040) (X(J,JNC),J=7,NVSB),X(NVAD,JNC) IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8040) X(NVAD,JNC) IF (JCST.EQ.0) WRITE(LUB,8040) (X(J,JNC),J=7,NVAD) 190 CONTINUE 200 IF(JPRT.EQ.0) GOTO 260 DO 220 JSPA=1,NVAR IF(JSPA.EQ.NVAR.AND.JCST.EQ.1) GOTO 220 JS=JSPA+1 DO 210 JSPB=JS,NVAD IF(JSPB.EQ.NVAR.AND.JCST.EQ.1) GOTO 210 JSM=JSPB-1 CALL CORR(X,JSPB,JSPA,NCAS,NVSB,JCST,AW,Y,RESDU,NDXX,NDXY, 1 SPEAR,PEARC,PREC) C(JSM,JSPA)=SPEAR H(JSM,JSPA)=PEARC 210 CONTINUE 220 CONTINUE WRITE(LUB,8080) LAB(NVAD) DIAG=1.0 WRITE(LUB,8090) LAB(1),DIAG DO 230 JSPA=1,NVSB JSPP=JSPA+1 IF(JCST.EQ.1.AND.JSPP.EQ.NVAR) GOTO 230 WRITE(LUB,8090) LAB(JSPP),(H(JSPA,JSPB),JSPB=1,JSPA),DIAG 230 CONTINUE IF(JCST.EQ.1) WRITE(LUB,8090) LAB(NVAD),(H(JSPA,JSPB), 1 JSPB=1,NVSB),DIAG IF (JCST.EQ.0) WRITE(LUB,8090) LAB(NVAD),(H(JSPA,JSPB), 1 JSPB=1,NVAR),DIAG DO 240 J=1,NVAR 240 JNDVC(J)=J WRITE(LUB,8100) LAB(NVAD) WRITE(LUB,8090) LAB(1),DIAG DO 250 JSPA=1,NVSB JSPP=JSPA+1 IF(JCST.EQ.1.AND.JSPP.EQ.NVAR) GOTO 250 WRITE(LUB,8090) LAB(JSPP),(C(JSPA,JSPB),JSPB=1,JSPA),DIAG 250 CONTINUE IF(JCST.EQ.1) WRITE(LUB,8090) LAB(NVAD),(C(JSPA,JSPB), 1 JSPB=1,NVSB),DIAG IF (JCST.EQ.0) WRITE(LUB,8090) LAB(NVAD),(C(JSPA,JSPB), 1 JSPB=1,NVAR),DIAG 260 DO 270 JNC=1,NCAS CC-----MEANWHILE, INITIALIZATION OF THE WEIGHTS WEIGHTS(JNC)=1.0 Y(JNC)=X(NVAD,JNC) 270 CONTINUE 8000 FORMAT(/' VARIABLE WITH LABEL ',A10,' IS CONSTANT OVER'/ 1 ' ALL OBSERVATIONS HENCE THE PROGRAM CANNOT RUN.', 1 /' PLEASE REPEAT THE ANALYSIS WITHOUT THIS VARIABLE.'/) 8010 FORMAT(//' MEDIANS = '/) 8015 FORMAT(//' DISPERSION OF ABSOLUTE VALUES = '/) 8020 FORMAT(7X,6(A10,1X)/15X,5(A10,1X)) 8030 FORMAT(7X,6(F10.4,1X)) 8040 FORMAT(15X,5(F10.4,1X)) 8050 FORMAT(//' DISPERSIONS = '/) 8060 FORMAT(//' THE STANDARDIZED OBSERVATIONS ARE:'/) 8070 FORMAT(I5,2X,6(F10.4,1X)) 8080 FORMAT(//' PEARSON CORRELATION COEFFICIENTS', 1 ' BETWEEN THE VARIABLES'/1X,'( ', 1 A10,' IS THE RESPONSE VARIABLE)'/) 8090 FORMAT(1X,A10,(2X,11F6.2)/) 8100 FORMAT(//' SPEARMAN RANK CORRELATION COEFFICIENTS', 1 ' BETWEEN THE VARIABLES'/1X,'( ',A10, 1 ' IS THE RESPONSE VARIABLE)'/ ) RETURN END CC CC *CORR* : CALCULATES THE SPEARMAN RANK CORRELATION COEFFICIENT CC BETWEEN VARIABLES JSPA AND JSPB. SUBROUTINE CORR(X,JSPB,JSPA,NCAS,NVSB,JCST,AW,Y,RESDU,NDXX,NDXY, 1 SPEAR,PEARC,PREC) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION X(NDXX,NDXY),AW(NDXY),Y(NDXY),RESDU(NDXY) INTEGER JCST,JSPA,JSPB,JTA,K,KK,NCM,NDX,NDXX,NDXY,NCAS INTEGER NJ,NJA,NJAP,NJB,NJP,NTWE,NVSB REAL AL,CVVAR,PEARC,PREC,SPEAR,SUMA,SUMB,STDA,STDB,TA,TB,W REAL X(NDXX,NDXY),AW(NDXY),Y(NDXY),RESDU(NDXY) CC-----IN THIS SUBROUTINE THE VECTOR "RESDU" IS USED AS A WORKING CC-----VECTOR FOR CALCULATING THE CORRELATION MATRICES SPEAR=0.0 SUMA=0.0 SUMB=0.0 STDA=0.0 STDB=0.0 CVVAR=0.0 AL=NCAS IF (((JSPB-JSPA).EQ.1).OR.(JCST.EQ.1.AND.JSPA.EQ.NVSB)) 1 GOTO 5 GOTO 105 5 NTWE=0 DO 10 NJ=1,NCAS AW(NJ)=X(JSPA,NJ) RESDU(NJ)=NJ 10 CONTINUE NCM=NCAS-1 15 DO 20 NJ=1,NCM NJP=NJ+1 NJA=NJ DO 30 NJB=NJP,NCAS IF(AW(NJB).GE.AW(NJA)) GOTO 30 NJA=NJB 30 CONTINUE IF(NJA.EQ.NJ) GOTO 20 W=AW(NJ) AW(NJ)=AW(NJA) AW(NJA)=W W=RESDU(NJA) RESDU(NJA)=RESDU(NJ) RESDU(NJ)=W 20 CONTINUE DO 40 NJ=1,NCAS NDX=RESDU(NJ) IF(NTWE.EQ.0) Y(NDX)=NJ IF(NTWE.NE.0) AW(NDX)=NJ 40 CONTINUE NJ=1 55 TA=1 NJA=RESDU(NJ) IF(NTWE.EQ.0) TB=Y(NJA) IF(NTWE.NE.0) TB=AW(NJA) 60 NJA=RESDU(NJ) NJAP=RESDU(NJ+1) IF(NTWE.EQ.0.AND.(X(JSPA,NJA).NE.X(JSPA,NJAP))) GOTO 70 IF(NTWE.NE.0.AND.(X(JSPB,NJA).NE.X(JSPB,NJAP))) GOTO 70 TA=TA+1.0 IF(NTWE.EQ.0) TB=TB+Y(NJAP) IF(NTWE.NE.0) TB=TB+AW(NJAP) NJ=NJ+1 IF (NJ.EQ.NCAS) GOTO 70 GOTO 60 70 IF(TA.EQ.1.0) GOTO 100 JTA=TA NJB=NJ+1-JTA DO 80 K=NJB,NJ KK=RESDU(K) IF(NTWE.EQ.0) Y(KK)=TB/TA IF(NTWE.NE.0) AW(KK)=TB/TA 80 CONTINUE 100 NJ=NJ+1 IF (NJ.LT.NCAS) GOTO 55 IF(NTWE.NE.0) GOTO 150 105 DO 110 NJ=1,NCAS AW(NJ)=X(JSPB,NJ) RESDU(NJ)=NJ SUMA=SUMA+X(JSPA,NJ) SUMB=SUMB+X(JSPB,NJ) 110 CONTINUE SUMA=SUMA/AL SUMB=SUMB/AL NTWE=1 GOTO 15 150 DO 120 NJ=1,NCAS SPEAR=SPEAR+(Y(NJ)-AW(NJ))**2 STDA=STDA+(X(JSPA,NJ)-SUMA)**2 STDB=STDB+(X(JSPB,NJ)-SUMB)**2 CVVAR=CVVAR+(X(JSPA,NJ)-SUMA)*(X(JSPB,NJ)-SUMB) 120 CONTINUE SPEAR=1.0-6.0*SPEAR/(AL*(AL**2-1.0)) STDA=SQRT(STDA*STDB) IF(ABS(STDA).GT.PREC) PEARC=CVVAR/(STDA) IF(ABS(STDA).LE.PREC) PEARC=99.99 RETURN END CC CC *PULL* : SEARCHES THE KTH VALUE IN A VECTOR CC OF LENGTH 'NCAS'. FUNCTION PULL(AW,NDXY,AA,NCAS,K) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION AA(NCAS),AW(NDXY) INTEGER J,JNC,K,L,LR,NCAS,NDXY REAL AX,PULL,WA REAL AA(NCAS),AW(NDXY) DO 10 JNC=1,NCAS 10 AW(JNC)=AA(JNC) L=1 LR=NCAS 20 IF(L.GE.LR) GOTO 90 AX=AW(K) JNC=L J=LR 30 IF(JNC.GT.J) GOTO 80 40 IF(AW(JNC).GE.AX) GOTO 50 JNC=JNC+1 GOTO 40 50 IF(AW(J).LE.AX) GOTO 60 J=J-1 GOTO 50 60 IF(JNC.GT.J) GOTO 70 WA=AW(JNC) AW(JNC)=AW(J) AW(J)=WA JNC=JNC+1 J=J-1 70 GOTO 30 80 IF(J.LT.K) L=JNC IF (K.LT.JNC) LR=J GOTO 20 90 PULL=AW(K) RETURN END CC CC *AMDAN* : CALCULATES THE MEDIAN OF A VECTOR CC OF LENGTH 'NCAS'. FUNCTION AMDAN(AW,NDXY,AA,NCAS) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION AA(NCAS),AW(NDXY) INTEGER JNDL,NCAS,NDXY REAL AL,AMDAN,PULL REAL AA(NCAS),AW(NDXY) AL=NCAS JNDL=INT(AL/2.0) IF (MOD(NCAS,2).EQ.0) THEN AMDAN= (PULL(AW,NDXY,AA,NCAS,JNDL)+ 1 PULL(AW,NDXY,AA,NCAS,JNDL+1))/2.0 ELSE AMDAN= PULL(AW,NDXY,AA,NCAS,JNDL+1) ENDIF RETURN END CC CC *RSQU* : CALCULATES COEFFICIENT OF DETERMINATION AND F-VALUE CC OF LS OR RLS CC FUNCTION RSQU(NCAS,NVAD,JCST,Y,NDXY,SSE,FVALUE,PREC, 1 XMAD,XMED,NDXX,WEIGHTS,NNUL) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION Y(NDXY),XMAD(NDXX),XMED(NDXX),WEIGHTS(NDXY) INTEGER JCST,JNC,NCAS,NDXX,NDXY,NNUL,NVAD REAL AL,DF1,DF2,FVALUE,PREC,RSQU,RYM,SSE,SSEL,SST REAL Y(NDXY),XMAD(NDXX),XMED(NDXX),WEIGHTS(NDXY) RYM=0.0 IF (JCST.EQ.0) GOTO 20 DO 10 JNC=1,NCAS 10 RYM=(Y(JNC)*XMAD(NVAD)+XMED(NVAD))*WEIGHTS(JNC)+RYM AL=NNUL RYM=RYM/AL 20 SST=0.0 DO 30 JNC=1,NCAS 30 SST=SST+(((Y(JNC)*XMAD(NVAD)+XMED(NVAD))-RYM)**2)*WEIGHTS(JNC) DF1=(NVAD-1)-JCST DF2=NNUL-(NVAD-1) IF(SST.LT.PREC)SST=PREC RSQU=1.0-(SSE/SST) IF (RSQU.LT.0.0) RSQU=0.0 IF (RSQU.GT.1.0) RSQU=1.0 SSEL=SSE IF(SSEL.LT.PREC)SSEL=PREC FVALUE=((SST-SSEL)/DF1)/(SSEL/DF2) IF(FVALUE.LT.(0.0))FVALUE=0.0 RETURN END CC CC *RESTT* : CALCULATES FOR ALL JNC (JNC=1,..,NCAS) THE RESIDUAL CC Y(JNC)-SUM (GESCF(J)*X(J,JNC)), FOR J=1,..,LTSTE. SUBROUTINE RESTT(GESCF,JABS,JTR,LTSTE,NCAS,NVAD,NZWE, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NDXX,NDXY,NOPT,SCAL,QUAN,PREC) C INCLUDE 'MID_REL_INCL:implicit.inc' C DIMENSION GESCF(LTSTE),X(NDXX,NDXY),Y(NDXY),WEIGHTS(NDXY) C DIMENSION RESDU(NDXY),XMED(NDXX),XMAD(NDXX) INTEGER J,JABS,JNC,JTR,LTSTE,NCAS,NDXX,NDXY,NOPT,NVAD,NZWE REAL ANVAR,AVLQS,PREC,QUAN,SCAL,WSC REAL GESCF(LTSTE),X(NDXX,NDXY),Y(NDXY),WEIGHTS(NDXY) REAL RESDU(NDXY),XMED(NDXX),XMAD(NDXX) CC-----JTR = 1, THEN CALCULATE THE RESIDUALS FOR THE STANDARDIZED OBSERVATIONS CC----- 0 " " " RAW CC-----JABS = 1, CALCULATE THE ABSOLUTE RESIDUALS CC----- 0 RESIDUALS CC-----NOPT = 1, CALCULATE THE LMS FINAL SCALE AND THE RESULTING WEIGHTS CC----- 0, OTHERWISE NZWE=0 SCAL=0.0 AVLQS=0.0 ANVAR=NVAD-1 DO 10 JNC=1,NCAS IF (JTR.NE.1) THEN RESDU(JNC)=Y(JNC) ELSE RESDU(JNC)=Y(JNC)*XMAD(NVAD)+XMED(NVAD) ENDIF DO 20 J=1,LTSTE IF (JTR.NE.1) THEN RESDU(JNC)=RESDU(JNC)-X(J,JNC)*GESCF(J) ELSE RESDU(JNC)=RESDU(JNC)-((X(J,JNC)*XMAD(J)+XMED(J))*GESCF(J)) ENDIF 20 CONTINUE IF (JABS.EQ.1.AND.NOPT.NE.1) RESDU(JNC)=ABS(RESDU(JNC)) IF (ABS(RESDU(JNC)).LT.PREC) NZWE=NZWE+1 IF (NOPT.EQ.0) GOTO 10 IF (ABS(RESDU(JNC)).LE.(2.5*QUAN)) THEN WSC=1.0 ELSE WSC=0.0 ENDIF AVLQS=AVLQS+WSC SCAL=SCAL+(RESDU(JNC)*RESDU(JNC))*WSC 10 CONTINUE C write(6,*) 'In PRC/RESTT: NOPT = ',NOPT IF (NOPT.EQ.0) RETURN SCAL=SQRT(SCAL/(AVLQS-ANVAR)) + PREC NZWE=0 DO 40 JNC=1,NCAS IF (ABS(RESDU(JNC)).LE.(2.5*SCAL)) THEN WEIGHTS(JNC)=1.0 NZWE=NZWE+1 ELSE WEIGHTS(JNC)=0.0 ENDIF C write(6,*) 'In PRC: ',jnc,resdu(jnc),scal,nzwe,weights(jnc) 40 CONTINUE RETURN END