C @(#)progress.for 17.1.1.1 (ES0-DMD) 01/25/02 17:56:14 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 PROGRESS : PROGRAM FOR ROBUST REGRESSION CC----------------------------------------------------------------------------- CC-----LIST OF VARIABLES AND ARRAYS IN THIS PROGRAM CC----------------------------------------------------------------------------- CC NDXY = MAXIMAL NUMBER OF CASES CC NMXV = MAXIMAL NUMBER OF EXPLANATORY VARIABLES CC NDXX = NMXV+1 = MAXIMAL NUMBER OF COEFFICIENTS, CC INCLUDING INTERCEPT CC JDMB = NMXV*(NMXV+1) CC----------------------------------------------------------------------------- CC X(NDXX,NDXY) = MATRIX OF THE VALUES OF THE EXPLANATORY CC VARIABLES (TRANSPOSED DESIGN MATRIX) CC Y(NDXY) = VECTOR WITH THE VALUES OF CC THE RESPONSE VARIABLE CC RESDU(NDXY) = VECTOR FOR STORING THE RESIDUALS CC (EXCEPT IN SUBROUTINE STATIS, WHERE IT IS CC A WORKING VECTOR) CC WEIGHTS(NDXY) = WEIGHT ASSIGNED TO EACH OBSERVATION CC AW(NDXY) = WORKING VECTOR CC A(NMXV+1) = VECTOR OF REGRESSION COEFFICIENTS (LS OR RLS) CC DA(NMXV+1) = VECTOR WITH THE STAND. DEVIATION OF THE COEFFICIENTS CC DRCTA(NMXV+1) = VECTOR OF COEFFICIENTS FOR THE LMS CC XMED(NMXV+1) = VECTOR WITH THE MEDIANS OF ALL VARIABLES CC XMAD(NMXV+1) = VECTOR WITH MEDIAN ABSOLUTE DEVIATIONS CC UHYB(NDXY) = VECTOR WITH THE RESISTANT DIAGNOSTICS CC C(NMXV,NMXV+1) = MATRIX USED FOR - STORING THE COEFFICIENTS CC OF THE LINEAR SYSTEM TO BE SOLVED CC AT EACH REPLICATION CC - STORING THE SPEARMAN CC CORRELATION COEFFICIENTS CC H(NMXV,NMXV+1) = MATRIX USED FOR - STORING THE VARIANCE-COVARIANCE CC MATRIX OF THE LS AND RLS COEFFICIENTS CC - STORING THE PEARSON CC CORRELATION COEFFICIENTS CC HVEC(JDMB) = AUXILIARY VECTOR CC----------------------------------------------------------------------------- CC NMVAL(NDXY) = VECTOR WITH THE ORIGINAL INDEX OF THE CASES CC JNDVC(NMXV+1) = VECTOR - CONTAINING THE INDEX OF THE CASES CC CONSIDERED IN A RANDOM SUBSAMPLE CC - COUNTING THE NUMBER OF MISSING CASES FOR CC EACH VARIABLE CC VALMS(NMXV+1) = VECTOR WITH THE "MISSING VALUE CODES" CC JNDEX(NMXV+1) = JNDEX(J), (J=1,NVAR+1) IS 1 IF VAR(J) CONTAINS MISSING CC VALUES AND IS 0 OTHERWISE CC (CONSIDERED IN A RANDOM SUBSAMPLE) CC JNDA(NMXV) = VECTOR CONTAINING THE INDEX OF THE CASES CC WHICH YIELD THE "BEST" LMS SOLUTION CC JPLACE(NMXV+1) = VECTOR CONTAINING THE POSITIONS OF THE VARIABLES CC USED IN THE ANALYSIS (FROM THE TOTAL DATA SET) CC LAB(NMXV+1) = VECTOR WITH THE LABEL OF EACH VARIABLE CC----------------------------------------------------------------------------- CC-----LIST OF MOST IMPORTANT PARAMETERS: CC----------------------------------------------------------------------------- CC FNAMEA = NAME OF THE FILE WHERE THE INPUT HAS TO BE READ CC FNAMEB = NAME OF THE FILE WHERE THE OUTPUT HAS TO BE WRITTEN CC FNAMEC = NAME OF THE FILE WHERE THE DATA HAS TO BE SAVED CC LUA = LOGICAL UNIT FOR 'FNAMEA' CC LUB = LOGICAL UNIT FOR 'FNAMEB' CC LUC = LOGICAL UNIT FOR 'FNAMEC' CC NVAR = NUMBER OF COEFFICIENTS (INCLUDING THE CONSTANT TERM CC IF THERE IS ONE) CC ANVAR = REAL VALUE OF NVAR CC NVSB = NVAR-1 CC NVAD = NVAR+1 CC NCAS = NUMBER OF CASES CC AL = REAL VALUE OF NCAS CC ALGO = CHARACTER WHICH IS 'E' FOR THE EXTENSIVE SEARCH CC VERSION AND 'Q' FOR THE QUICK VERSION OF THE ALGORITHM CC NZWE = NUMBER OF NON-ZERO WEIGHTS CC RSQ = R SQUARED (COEFFICIENT OF DETERMINATION) CC AVW = AVERAGE WEIGHT CC FCKW = OBJECTIVE FUNCTION VALUE CC JCST = INTEGER WHICH IS 1 IN THE CASE OF A MODEL WITH CC INTERCEPT AND 0 OTHERWISE CC JPLT = PLOT OPTION (0=NONE, 1=RESIDUAL, 2=INDEX, 3=BOTH) CC JPRT = PRINT OPTION (0=SMALL, 1=MEDIUM-SIZED, 2=LARGE) CC MVAL = OPTION FOR MISSING VALUES (0=NONE, 1=DELETE, 2=REPLACE) CC YNSAVE = CHARACTER (YES OR NO) NEEDED FOR THE INTERACTIVE INPUT CC JHEAD = TITLE FOR THE OUTPUT CC NZWE = NUMBER OF NON-ZERO WEIGHTS CC PREC = PRECISION OF THE CALCULATIONS CC JEY,JLP = INTEGER VARIABLES USED FOR COUNTING THE NUMBER OF CC REPLICATIONS IN THE LMS ALGORITHM CC JPLXY = OPTION FOR XY PLOT (0=NO XY PLOT, 1=OTHERWISE) CC JREP = NUMBER OF RANDOM SUBSAMPLES TO BE CONSIDERED CC JREG = WHICH REGRESSION ? (1=LS, 2=LMS, 3=RLS) CC JDIAG = OPTION FOR DIAGNOSTICS (0=NO DIAGNOSTICS, ELSE 1) CC JQU = THE RANK OF THE ABSOLUTE RESIDUAL TO BE MINIMIZED CC NSTOP = 1 MEANS STOP THE ANALYSIS, OTHERWISE NSTOP=0 CC----------------------------------------------------------------------------- CC-----LIST OF SUBROUTINES AND FUNCTIONS CC----------------------------------------------------------------------------- CC SUBROUTINE WRDATA :CREATES FIRST PART OF THE OUTPUT CC SUBROUTINE WRDATB :CREATES SECOND PART OF THE OUTPUT CC SUBROUTINE SUBHAT :CALCULATION OF THE HAT MATRIX CC SUBROUTINE PRTRLS :PRINTS THE RLS (OR LS) RESULTS CC SUBROUTINE SUBLMS :ALGORITHM FOR CALCULATING THE LMS ESTIMATES CC SUBROUTINE PRTLMS :PRINTS THE LMS RESULTS CC SUBROUTINE RDAT :INTERACTIVE INPUT WITH READING OF THE DATA CC FUNCTION JNTGR :TRANSFORMS A CHARACTER INTO AN INTEGER CC SUBROUTINE SMISSING:HANDLES MISSING VALUES CC SUBROUTINE SUBREP :CALCULATES THE NUMBER OF REPLICATIONS CC FOR THE LMS ALGORITHM CC SUBROUTINE STATIS :STANDARDIZATION OF THE OBSERVATIONS CC AND CALCULATION OF SOME STATISTICS CC SUBROUTINE CORR :CALCULATES SPEARMAN RANK CORRELATION COEFFICIENT CC FUNCTION PULL :SEARCHES THE KTH VALUE IN A VECTOR CC FUNCTION AMDAN :CALCULATES THE MEDIAN CC FUNCTION RSQU :CALCULATES COEFFICIENT OF DETERMINATION CC OF LS AND RLS CC SUBROUTINE RESTT :CALCULATES RESIDUALS AND LMS SCALE CC SUBROUTINE RANGS :SORTS A VECTOR CC SUBROUTINE RDUAL :GENERATES A TABLE WITH OBSERVED AND ESTIMATED CC RESPONSE, RESIDUALS AND STANDARDIZED RESIDUALS CC AND (IF WANTED) CALLS RESIDUAL PLOTS CC SUBROUTINE SHHLF :CALCULATES LMS IN ONE DIMENSION CC SUBROUTINE TRC :RETRANSFORMATION OF THE VARIANCE-COVARIANCE MATRIX CC SUBROUTINE SCHCV :PRINTS THE VARIANCE-COVARIANCE MATRIX CC SUBROUTINE RTRAN :RETRANSFORMATION OF THE REGRESSION COEFFICIENTS CC SUBROUTINE RANPN :RANDOM NUMBER GENERATOR CC SUBROUTINE GENPN :GENERATES ALL COMBINATIONS OF P POINTS OUT OF N CC SUBROUTINE EQUAT :SOLVES A SYSTEM OF LINEAR EQUATIONS CC SUBROUTINE LSL :LS OR RLS IN THE CASE P=1 CC SUBROUTINE FCN :PUTS A ROW OF THE MATRIX X IN A VECTOR CC SUBROUTINE LSREG :CALCULATES LS OR RLS REGRESSION CC FUNCTION QLSRG :CALCULATES LS OR RLS OBJECTIVE FUNCTION CC SUBROUTINE MATNV :MATRIX INVERSION CC SUBROUTINE LCAT :TREATS THE LOCATION AND SCALE ESTIMATORS CC SUBROUTINE SMISLOC :TREATMENT OF MISSING VALUES IN THE CASE OF LOCATION CC SUBROUTINE MOVE :MOVES CHARACTERS TO THE RIGHT (IN A LABEL) CC FUNCTION PVAL :CALCULATES THE P-VALUE ASSOCIATED WITH AN F-VALUE CC FUNCTION PTVAL :CALCULATES THE P-VALUE ASSOCIATED WITH A T-VALUE CC FUNCTION SUM :FUNCTION NEEDED IN PVAL CC FUNCTION ODD :FUNCTION NEEDED IN PVAL CC SUBROUTINE GRAF :MAKES A SCATTERGRAM, AND RESIDUAL AND/OR INDEX PLOTS CC----------------------------------------------------------------------------- CC PROGRAM PROGRESS INCLUDE 'MID_REL_INCL:implicit.inc' CHARACTER YNSAVE,ALGO CHARACTER*30 FNAMEA,FNAMEB,FNAMEC CHARACTER*60 JHEAD CC ------------------------------------------------------------------------ CC IF YOU LIKE TO CHANGE THE DIMENSIONS OF THE ARRAYS, THEN YOU CC ONLY HAVE TO ADAPT THE VALUES IN THE FOLLOWING NINE LINES: CC ------------------------------------------------------------------------ DIMENSION X(21,3000),Y(3000),RESDU(3000) DIMENSION AW(3000),UHYB(3000) DIMENSION JNDVC(21),VALMS(21),JNDEX(21) DIMENSION NMVAL(3000),WEIGHTS(3000) DIMENSION A(21),DA(21),JNDA(20),JPLACE(21) DOUBLE PRECISION H(20,21),HVEC(420) DOUBLE PRECISION DTERM DIMENSION DRCTA(21),XMED(21),XMAD(21) DIMENSION C(20,21) CHARACTER*10 LAB(21) INTEGER OUTN1,OUTN2 INTEGER ISTAT CHARACTER*80 OUTEXT INCLUDE 'MID_INCLUDE:st_def.inc' COMMON /VMR/ MADRID INCLUDE 'MID_INCLUDE:st_dat.inc' NMXV=20 NDXY=3000 CC ------------------------------------------------------------------------ CC INITIALIZATIONS CC ------------------------------------------------------------------------ CALL STSPRO('progress') LUA=1 LUB=2 LUC=3 NDXX=NMXV+1 JDMB=NMXV*NDXX NZWE=0 PREC=1.0E-6 JLP=0 JPLXY=0 JREG=0 NSTOP=0 DO 5 J=1,NMXV JNDA(J)=0 5 CONTINUE CC ------------------------------------------------------------------------ CC READ THE DATA CC ------------------------------------------------------------------------ CALL RMIDAT(NCAS,NVAR,JCST,JPRT,NVSB,NVAD,VALMS, 1 X,Y,JPLT,JNDEX,PREC,JPLACE,LAB,XMED,XMAD,AW,NMVAL, 1 RESDU,WEIGHTS,NDXX,NDXY,MVAL,NMXV,LUA,LUB,LUC,JHEAD, 1 YNSAVE,JDIAG,ALGO,NSTOP,FNAMEA,FNAMEB,FNAMEC,OUTN1,OUTN2) IF (NSTOP.EQ.1) GOTO 30 AL=NCAS ANVAR=NVAR CC ------------------------------------------------------------------------ CC CREATE THE BEGINNING OF THE OUTPUT: THE SELECTED OPTIONS CC ------------------------------------------------------------------------ CALL WRDATA(NVAR,NCAS,JCST,MVAL,JPRT,ALGO,YNSAVE,FNAMEA,FNAMEB, 1 FNAMEC,JHEAD,LUB) CC ------------------------------------------------------------------------ CC IF THE DATA CONTAIN MISSING VALUES, THEN TREAT THEM CC ------------------------------------------------------------------------ IF (MVAL.NE.0) THEN CALL SMISSING(NVAR,NCAS,NDXX,NDXY,JCST,X,Y,AW,NMVAL,MVAL, 1 JNDEX,JNDVC,VALMS,NSTOP,LAB,JPRT,LUB) IF (NSTOP.EQ.1) GOTO 20 ENDIF CC ------------------------------------------------------------------------ CC CREATE SECOND PART OF THE OUTPUT: THE DATA AND AN X-Y PLOT CC ------------------------------------------------------------------------ CALL WRDATB(NVAR,NCAS,NDXX,NDXY,JCST,MVAL,LAB,NMVAL,X,Y, 1 AW,JPLT,JPRT,JREG,JHEAD,LUB) CC ------------------------------------------------------------------------ CC STANDARDIZATION OF THE DATA (WITH OPTIONAL PRINTING) CC AND CALCULATION OF SOME STATISTICS ON THE VARIABLES CC (IN THIS SUBROUTINE "RESDU" IS USED AS A WORKING VECTOR) CC ------------------------------------------------------------------------ CALL STATIS(X,Y,XMED,XMAD,RESDU,WEIGHTS,AW,NMVAL,C,H,JNDVC, 1 NSTOP,NVAR,NVAD,NVSB,NCAS,JCST,JPRT,NDXX,NDXY,LUB,NMXV, 1 PREC,LAB) IF (NSTOP.EQ.1) GOTO 20 CC ------------------------------------------------------------------------ CC CALCULATION OF THE LS ESTIMATES CC ------------------------------------------------------------------------ IF (NCAS.LE.(NVAR*2)) THEN CALL STTPUT('TOO MANY COEFFICIENTS ACCORDING TO THE ',ISTAT) CALL STTPUT('NUMBER OF CASES. NUMBER OF CASES',ISTAT) WRITE(OUTEXT,8050) NCAS,NVAR CALL STTPUT(OUTEXT,ISTAT) CALL STTPUT('THE NUMBER OF CASES MUST BE TWICE THE NUMBER 1 OF COEFFICIENTS.',ISTAT) IF (JCST.EQ.1) CALL STTPUT(' (INCLUDING THE 1 CONSTANT TERM!)',ISTAT) GOTO 20 ENDIF JREG=1 IF (NVAR.EQ.1) THEN CALL LSL(NCAS,NDXY,X,Y,WEIGHTS,A,FCKW,H,NMXV,NDXX) ELSE CALL LSREG(NDXX,NDXY,NMXV,NVAR,NCAS,NVAR,A,X,Y,WEIGHTS,DA, 1 H,FCKW,HVEC,JDMB,JNDEX) ENDIF CC ------------------------------------------------------------------------ CC RETRANSFORMATION OF THE LS ESTIMATES CC ------------------------------------------------------------------------ CALL RTRAN(NVAR,JCST,NVSB,NVAD,NDXX,XMED,XMAD,A,NVAD,FCKW) CALL TRC(H,DA,NMXV,NDXX,NDXY,NVAR,JCST,NVSB,NVAD,XMED,XMAD) CC ------------------------------------------------------------------------ CC PRINTING THE LS RESULTS WITH THE REQUESTED OPTIONS CC ------------------------------------------------------------------------ CALL PRTRLS(NVAR,NCAS,JCST,NVSB,NVAD,NDXX,NDXY,NMXV, 1 XMED,XMAD,A,DA,H,AW,RESDU,WEIGHTS,X,Y,NMVAL,AVW,NZWE,FCKW, 1 JREG,JPRT,JPLT,JDIAG,JHEAD,LAB,PREC,LUB) CC ------------------------------------------------------------------------ CC CALCULATION OF LMS REGRESSION CC ------------------------------------------------------------------------ C CALL STTPUT('',ISTAT) CALL STTPUT(' ALGORITHM FOR CALCULATING THE LEAST MEDIAN 1 OF SQUARES IS STARTED. ',ISTAT) CALL SUBLMS(NVAR,NCAS,JCST,NDXX,NDXY,NMXV,PREC,X,Y, 1 XMED,XMAD,JDIAG,ALGO,FNAMEB,JQU,JREP,DRCTA,JNDA,RESDU,WEIGHTS, 1 UHYB,AVW,NZWE,RSQ,JHYB,JDMB,JERD,JLP,NERR,DTERM,C,HVEC,JNDVC,AW) CC ------------------------------------------------------------------------ CC RETRANSFORMATION OF THE LMS ESTIMATES CC ------------------------------------------------------------------------ CALL RTRAN(NVAR,JCST,NVSB,NVAD,NDXX,XMED,XMAD,DRCTA,NVAD,FCKW) DRCTA(NVAD)=DRCTA(NVAD)*XMAD(NVAD) CC ------------------------------------------------------------------------ CC PRINTING THE LMS RESULTS WITH THE REQUESTED OPTIONS CC ------------------------------------------------------------------------ CALL PRTLMS(NVAR,NCAS,JCST,NVSB,NVAD,NDXX,NDXY,NMXV,JNDVC,XMED, 1 XMAD,DRCTA,AW,RESDU,WEIGHTS,X,Y,NMVAL,UHYB,AVW,NZWE,RSQ,JREG, 1 JPRT,JPLT,JREP,JERD,JLP,JQU,JHYB,JDIAG,NSTOP,JHEAD,LAB,PREC,LUB) IF (NSTOP.NE.1) THEN CC ------------------------------------------------------------------------ CC CALCULATION OF RLS REGRESSION CC ------------------------------------------------------------------------ JREG=3 CALL STTPUT('PROGRESS IS CURRENTLY CALCULATING REWEIGHTED LEAST 1 SQUARES ...',ISTAT) NNUL=NZWE DO 10 J=1,NDXX A(J)=0.0 DA(J)=0.0 DO 10 JNC=1,NMXV 10 H(JNC,J)=0.D0 IF (NVAR.EQ.1) THEN CALL LSL(NCAS,NDXY,X,Y,WEIGHTS,A,FCKW,H,NMXV,NDXX) ELSE CALL LSREG(NDXX,NDXY,NMXV,NVAR,NCAS,NVAR,A,X,Y,WEIGHTS,DA, 1 H,FCKW,HVEC,JDMB,JNDEX) ENDIF CC ------------------------------------------------------------------------ CC RETRANSFORMATION OF THE RLS ESTIMATES CC ------------------------------------------------------------------------ CALL RTRAN(NVAR,JCST,NVSB,NVAD,NDXX,XMED,XMAD,A,NVAD,FCKW) CALL TRC(H,DA,NMXV,NDXX,NDXY,NVAR,JCST,NVSB,NVAD,XMED,XMAD) CC ------------------------------------------------------------------------ CC PRINTING THE RLS RESULTS WITH THE REQUESTED OPTIONS CC ------------------------------------------------------------------------ CALL PRTRLS(NVAR,NCAS,JCST,NVSB,NVAD,NDXX,NDXY,NMXV, 1 XMED,XMAD,A,DA,H,AW,RESDU,WEIGHTS,X,Y,NMVAL,AVW,NZWE,FCKW, 1 JREG,JPRT,JPLT,JDIAG,JHEAD,LAB,PREC,LUB) CC-----UPDATES MIDAS TABLES CALL MIDUAL(A,JKEUS,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL,NDXX,NDXY, 1 JHEAD,LAB,LUA,OUTN1,OUTN2) ENDIF CC END OF THE RUN CALL STTPUT(' ',ISTAT) CALL STTPUT(' THE RUN HAS SUCCESSFULLY BEEN EXECUTED. ',ISTAT) CALL STTPUT(' ',ISTAT) 20 IF (YNSAVE.EQ.'Y') THEN WRITE (OUTEXT,8020) FNAMEC CALL STTPUT(OUTEXT,ISTAT) ENDIF IF (FNAMEA.NE.'CON') THEN WRITE(OUTEXT,8030) FNAMEA CALL STTPUT(OUTEXT,ISTAT) ENDIF IF (FNAMEB.NE.'CON'.AND.FNAMEB.NE.'PRN') THEN WRITE(OUTEXT,8040)FNAMEB CALL STTPUT(OUTEXT,ISTAT) ENDIF CC ----------------------------------------------------------------- CC FORMATS CC ----------------------------------------------------------------- 8020 FORMAT(' THE DATA IS SAVED IN FILE : ',A30) 8030 FORMAT(' THE DATA HAS BEEN READ FROM FILE : ',A30) 8040 FORMAT(' THE OUTPUT HAS BEEN WRITTEN IN FILE : ',A30) 8050 FORMAT( 8X,'= ',I5,' NUMBER OF COEFFICIENTS = ',I5) 30 CONTINUE CALL STSEPI END CC ------------------------------------------------------------------------ CC SUBROUTINE FOR CREATING THE BEGINNING OF THE OUTPUT: CC - DESCRIPTION OF THE SELECTED OPTIONS CC ------------------------------------------------------------------------ SUBROUTINE WRDATA(NVAR,NCAS,JCST,MVAL,JPRT,ALGO,YNSAVE,FNAMEA, 1 FNAMEB,FNAMEC,JHEAD,LUB) INCLUDE 'MID_REL_INCL:implicit.inc' CHARACTER ALGO,YNSAVE CHARACTER*30 FNAMEA,FNAMEB,FNAMEC CHARACTER*60 JHEAD WRITE(LUB,8050) WRITE(LUB,8040)JHEAD IF (JCST.EQ.0) THEN WRITE(LUB,8000) NCAS,NVAR ELSEIF (JCST.NE.0) THEN WRITE(LUB,8010) NCAS,NVAR ENDIF IF (ALGO.EQ.'Q') THEN WRITE(LUB,8020) ELSEIF (ALGO.EQ.'E') THEN WRITE(LUB,8030) ENDIF IF (FNAMEB.EQ.'CON') PAUSE ' ' CC-----TREATMENT OF MISSING VALUES IF (MVAL.EQ.0) THEN WRITE(LUB,8060) ELSEIF (MVAL.EQ.1) THEN WRITE(LUB,8070) ELSEIF (MVAL.EQ.2) THEN WRITE(LUB,8080) ENDIF IF (JPRT.EQ.0) THEN WRITE(LUB,8090) ELSEIF (JPRT.EQ.1) THEN WRITE(LUB,8100) ELSEIF (JPRT.EQ.2) THEN WRITE(LUB,8110) ENDIF IF (YNSAVE.EQ.'Y') THEN WRITE(LUB,8120) FNAMEC ENDIF IF (FNAMEA.NE.'CON') THEN WRITE(LUB,8130) FNAMEA ENDIF 8000 FORMAT(//' REGRESSION WITHOUT A CONSTANT TERM.'// 1 ' NUMBER OF CASES',8X,'= ',I5/ 1 ' NUMBER OF COEFFICIENTS = ',I5) 8010 FORMAT(/' REGRESSION WITH A CONSTANT TERM.'// 1 ' NUMBER OF CASES',8X,'= ',I5/ 1 ' NUMBER OF COEFFICIENTS (INCLUDING CONSTANT TERM) = ',I5) 8020 FORMAT(/' THE QUICK VERSION WILL BE USED.') 8030 FORMAT(/' THE EXTENSIVE SEARCH VERSION WILL BE USED.') 8040 FORMAT(/' DATA SET = ',A60) 8050 FORMAT(//30X,19('*')/30X,'* P R O G R E S S *'/ 1 30X,19('*')///' This program performs a robust regression', 1 ' analysis based on'/' the least median of squares', 1 ' (LMS) method as described in'/3X, 1 ' P. Rousseeuw (1984), Least Median of Squares Regression',/ 1 ' Journal of the American Statistical Association, 79,', 1' 871-880.'/' A user manual to this program is the book:'/ 1 ' P. Rousseeuw and A. Leroy (1987), Robust Regression'/ 1 ' and Outlier Detection, Wiley, New York.'/) 8060 FORMAT(/' THERE ARE NO MISSING VALUES.') 8070 FORMAT(/' TREATMENT OF MISSING VALUES IN OPTION 1:', 1 ' THIS MEANS THAT A CASE WITH A/14H MISSING VALUE', 1 ' FOR AT LEAST ONE VARIABLE WILL BE DELETED.'/) 8080 FORMAT(/' TREATMENT OF MISSING VALUES IN OPTION 2:'/ 1 ' FIRST, A CASE WITH A MISSING VALUE FOR THE RESPONSE VARIABLE' 1 /' OR FOR ALL EXPLANATORY VARIABLES WILL BE DELETED.'/ 1 ' THEN, A MISSING VALUE FOR A VARIABLE WILL BE REPLACED BY'/ 1 ' THE MEDIAN OF THE NON-MISSING VALUES .'/) 8090 FORMAT(/' SMALL OUTPUT IS WANTED.') 8100 FORMAT(/' MEDIUM-SIZED OUTPUT IS WANTED.') 8110 FORMAT(/' LARGE OUTPUT IS WANTED.') 8120 FORMAT(//' THE DATA ARE SAVED IN FILE : ',A30/) 8130 FORMAT(/' YOUR DATA RESIDE IN FILE',7X,': ',A30/) RETURN END CC ------------------------------------------------------------------------ CC SUBROUTINE FOR CREATING THE SECOND PART OF THE OUTPUT: CC - WRITING DOWN THE OBSERVATIONS CC - AN X-Y PLOT CC ------------------------------------------------------------------------ SUBROUTINE WRDATB(NVAR,NCAS,NDXX,NDXY,JCST,MVAL,LAB,NMVAL,X,Y, 1 AW,JPLT,JPRT,JREG,JHEAD,LUB) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION X(NDXX,NDXY),Y(NDXY),AW(NDXY),NMVAL(NDXY) CHARACTER*10 LAB(NDXX) CHARACTER*60 JHEAD NVSB=NVAR-1 NVAD=NVAR+1 IF ((NVAR.EQ.1.AND.JCST.EQ.0).OR. 1 (NVAR.EQ.2.AND.JCST.EQ.1)) JPLXY=1 IF (JPRT.EQ.2) THEN IF (MVAL.EQ.0) WRITE(LUB,8110) IF (MVAL.NE.0) WRITE(LUB,8120) IF (JCST.EQ.1) WRITE(LUB,8130) (LAB(J),J=1,NVSB),LAB(NVAD) IF (JCST.EQ.0) WRITE(LUB,8130) (LAB(J),J=1,NVAD) ENDIF DO 10 JNC=1,NCAS IF(JPRT.GE.2) THEN IF((NVAD.GT.6.AND.JCST.EQ.0).OR.(NVAD.GT.7.AND.JCST.EQ.1)) 1 GOTO 15 IF (JCST.EQ.1) WRITE(LUB,8140) NMVAL(JNC),(X(J,JNC),J=1,NVSB), 1 X(NVAD,JNC) IF (JCST.EQ.0) WRITE(LUB,8140) NMVAL(JNC),(X(J,JNC),J=1,NVAD) GOTO 5 15 WRITE(LUB,8140) NMVAL(JNC),(X(J,JNC),J=1,6) IF (JCST.EQ.1.AND.NVSB.GE.7) 1 WRITE(LUB,8150) (X(J,JNC),J=7,NVSB),X(NVAD,JNC) IF (JCST.EQ.1.AND.NVSB.LT.7) WRITE(LUB,8150) X(NVAD,JNC) IF (JCST.EQ.0) WRITE(LUB,8150) (X(J,JNC),J=7,NVAD) ENDIF 5 IF (JPLXY.EQ.1) THEN AW(JNC)=X(1,JNC) Y(JNC)=X(NVAD,JNC) ENDIF 10 CONTINUE IF (JPLXY.EQ.1) 1 CALL GRAF(AW,Y,NCAS,JPLXY,0,JPLT,NMVAL,NDXY,LUB,JREG, 1 JHEAD,NDXX,NVAD,LAB) 8110 FORMAT(/' THE OBSERVATIONS ARE:'/) 8120 FORMAT(/' THE OBSERVATIONS, AFTER TREATMENT' 1 ' OF MISSING VALUES ARE:'/) 8130 FORMAT(7X,6(A10,1X)/15X,5(A10,1X)) 8140 FORMAT(I5,2X,6(F10.4,1X)) 8150 FORMAT(15X,5(F10.4,1X)) RETURN END CC ------------------------------------------------------------------------ CC SUBROUTINE FOR COMPUTING THE ELEMENTS OF THE HAT MATRIX CC ------------------------------------------------------------------------ SUBROUTINE SUBHAT(NVAR,NCAS,NDXX,NDXY,NMXV,X,XMED,XMAD,H,AW, 1 WEIGHTS,NMVAL,FCKW,LUB) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION X(NDXX,NDXY),AW(NDXY),NMVAL(NDXY),WEIGHTS(NDXY) DOUBLE PRECISION H(NMXV,NDXX) DIMENSION XMED(NDXX),XMAD(NDXX) DO 10 J=1,NVAR DO 10 J2=1,J H(J,J2)=H(J,J2)/DBLE(FCKW) H(J2,J)=H(J,J2) 10 CONTINUE IF (NCAS.LE.20) WRITE (LUB,8000) IF (NCAS.GT.20) WRITE (LUB,8010) DO 20 JNA=1,NCAS DO 30 JNB=1,NCAS IF (NCAS.GT.20.AND.JNA.NE.JNB) GOTO 30 AW(JNB)=0.0 DO 40 J=1,NVAR DO 50 J2=1,NVAR RHJ2J=H(J2,J) AW(JNB)=AW(JNB)+(X(J2,JNA)*XMAD(J2)+XMED(J2))*WEIGHTS(JNA)* 1 RHJ2J*(X(J,JNB)*XMAD(J)+XMED(J))*WEIGHTS(JNB) 50 CONTINUE 40 CONTINUE 30 CONTINUE IF (NCAS.LE.9) THEN WRITE(LUB,8020) NMVAL(JNA),(AW(JNB),JNB=1,NCAS) GOTO 20 ENDIF IF (NCAS.GT.20) THEN WRITE(LUB,8020) NMVAL(JNA),AW(JNA) ELSE WRITE(LUB,8020) NMVAL(JNA),(AW(JNB),JNB=1,9) WRITE(LUB,8030) (AW(JNB),JNB=10,NCAS) ENDIF 20 CONTINUE 8000 FORMAT(//' --- THE HAT MATRIX ---'/1X,22('-')) 8010 FORMAT(//' THE DIAGONAL ELEMENTS OF THE HAT MATRIX'/1X,40('-')) 8020 FORMAT(1X,I4,1X,9(F6.3,2X)) 8030 FORMAT(4X,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3,F8.3) RETURN END CC ------------------------------------------------------------------------ CC PRINTING THE RLS ESTIMATES AND SOME OPTIONS CC ------------------------------------------------------------------------ SUBROUTINE PRTRLS(NVAR,NCAS,JCST,NVSB,NVAD,NDXX,NDXY,NMXV, 1 XMED,XMAD,A,DA,H,AW,RESDU,WEIGHTS,X,Y,NMVAL,AVW,NZWE,FCKW, 1 JREG,JPRT,JPLT,JDIAG,JHEAD,LAB,PREC,LUB) INCLUDE 'MID_REL_INCL:implicit.inc' INTEGER UNIT,STATUS DIMENSION X(NDXX,NDXY),Y(NDXY),WEIGHTS(NDXY),A(NMXV),DA(NMXV) DIMENSION XMED(NDXX),XMAD(NDXX),AW(NDXY),RESDU(NDXY),NMVAL(NDXY) DOUBLE PRECISION H(NMXV,NDXX),DS,PTVAL,DST,PVAL CHARACTER*10 LAB(NDXX) CHARACTER*60 JHEAD AL=NCAS ANVAR=NVAR WRITE(LUB,8000) IF (JREG.EQ.1) THEN WRITE(LUB,8010) NNUL=NCAS AVW=1 ELSE WRITE(LUB,8020) NNUL=NZWE ENDIF NDF=NNUL-NVAR CC-----PRINTING THE COEFFICIENTS WITH THEIR STANDARD ERROR, CC-----AND T-VALUE WRITE(LUB,8030) C DO 10 J=1,NVSB IF (DA(J).LT.PREC) DA(J)=PREC STUD=A(J)/DA(J) DS=DBLE(STUD) DST=PTVAL(DS,NDF) RDST=DST WRITE(LUB,8040) LAB(J),A(J),DA(J),STUD,RDST C----put value into MIDAS keyword OUTPUTR(1): CALL STKWRR('OUTPUTR',A(J),J,1,UNIT,STATUS) C----put also error into MIDAS keyword, use OUTPUTR(3) CALL STKWRR('OUTPUTR',DA(J),J+2,1,UNIT,STATUS) C 10 CONTINUE IF(DA(NVAR).LT.PREC) DA(NVAR)=PREC STUD=A(NVAR)/DA(NVAR) DS=DBLE(STUD) DST=PTVAL(DS,NDF) RDST=DST IF (JCST.EQ.0) WRITE(LUB,8040) LAB(NVAR),A(NVAR),DA(NVAR), 1 STUD,RDST IF (JCST.EQ.1) WRITE(LUB,8050) A(NVAR),DA(NVAR),STUD,RDST C----put value into MIDAS keyword OUTPUTR(2): CALL STKWRR('OUTPUTR',A(NVAR),NVAR,1,UNIT,STATUS) C----put error into MIDAS keyword, use OUTPUTR(4) CALL STKWRR('OUTPUTR',DA(NVAR),NVAR+2,1,UNIT,STATUS) C CC-----CALCULATION OF THE COEFFICIENT OF DETERMINATION AND F-VALUE RSQ=RSQU(NCAS,NVAD,JCST,Y,NDXY,FCKW,FVALUE,PREC, 1 XMAD,XMED,NDXX,WEIGHTS,NNUL) CC-----CALCULATION OF THE SCALE ESTIMATE A(NVAD) = SQRT(FCKW/(AL*AVW-ANVAR)) IF (JREG.EQ.1) THEN WRITE(LUB,8080) FCKW ELSE WRITE(LUB,8090) FCKW ENDIF WRITE(LUB,8060) NDF WRITE(LUB,8070) A(NVAD) CC-----PRINTING THE VARIANCE-COVARIANCE MATRIX CALL SCHCV(NVAR,NDXX,NMXV,H,LUB) CC-----IF DIAGNOSTICS ARE REQUESTED, THEN PRINT THE HAT MATRIX CC----- OR THE DIAGONAL ELEMENTS ONLY FCKW=FCKW/(AL*AVW-ANVAR) IF(FCKW.LT.PREC)GOTO 24 IF (JDIAG.EQ.1) CALL SUBHAT(NVAR,NCAS,NDXX,NDXY,NMXV,X,XMED, 1 XMAD,H,AW,WEIGHTS,NMVAL,FCKW,LUB) CC-----PRINT COEFFICIENT OF DETERMINATION, F-VALUE, AND P-VALUE 24 WRITE(LUB,8100) RSQ NDF1=NVAR-JCST DS=DBLE(FVALUE) DST=PVAL(DS,NDF1,NDF) RDST=DST WRITE(LUB,8110) FVALUE,NDF1,NDF,RDST IF (JREG.EQ.3) THEN WRITE(LUB,8200) NNUL WRITE(LUB,8210) AVW JKEUS=2 ELSE JKEUS=0 ENDIF CC-----GIVE RESIDUAL PLOTS IF REQUESTED IF(JPRT.NE.0.OR.JPLT.NE.0) THEN CALL RDUAL(A,JKEUS,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL,NDXX,NDXY, 1 JHEAD,LAB) ENDIF 8000 FORMAT(/1X,78('*')/) 8010 FORMAT(' LEAST SQUARES REGRESSION '/1X,24('*')/) 8020 FORMAT(' REWEIGHTED LEAST SQUARES BASED ON THE LMS'/ 1 1X,41('*')/) 8030 FORMAT(/5X,'VARIABLE',5X,'COEFFICIENT', 1 4X,'STAND. ERROR',5X,'T - VALUE',5X,'P - VALUE'/ 1 3X,70('-')) 8040 FORMAT(3X,A10,5X,F11.5,4X,F12.5,5X,F9.5,5X,F9.5) 8050 FORMAT(5X,'CONSTANT',5X,F11.5,4X,F12.5,5X,F9.5,5X,F9.5) 8060 FORMAT(/' DEGREES OF FREEDOM ',5X,'= ',4X,I5) 8070 FORMAT(/' SCALE ESTIMATE',10X,'= ',F15.5/) 8080 FORMAT(/' SUM OF SQUARES',10X,'= ',F15.5) 8090 FORMAT(/' WEIGHTED SUM OF SQUARES'' = ',F15.5) 8100 FORMAT(/' COEFFICIENT OF DETERMINATION (R SQUARED) =',F15.5/) 8110 FORMAT(' THE F-VALUE = ',F12.3,' (WITH ',I3,' AND ', 1 I4,' DF)',3X,'P - VALUE = ',F7.5) 8200 FORMAT(/'THERE ARE',2X,I4, 1 ' POINTS WITH NON-ZERO WEIGHT.' ) 8210 FORMAT(/' AVERAGE WEIGHT',10X,'= ',F15.5/) RETURN END CC ------------------------------------------------------------------------ CC SUBROUTINE FOR CALCULATING THE LMS REGRESSION CC ------------------------------------------------------------------------ SUBROUTINE SUBLMS(NVAR,NCAS,JCST,NDXX,NDXY,NMXV,PREC,X,Y, 1 XMED,XMAD,JDIAG,ALGO,FNAMEB,JQU,JREP,DRCTA,JNDA,RESDU,WEIGHTS, 1 UHYB,AVW,NZWE,RSQ,JHYB,JDMB,JERD,JLP,NERR,DTERM,C,HVEC,JNDVC,AW) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION X(NDXX,NDXY),Y(NDXY),C(NMXV,NDXX),DRCTA(NDXX) DIMENSION JNDA(NMXV),JNDVC(NDXX),XMED(NDXX),XMAD(NDXX) DIMENSION AW(NDXY),RESDU(NDXY),UHYB(NDXY),WEIGHTS(NDXY) DOUBLE PRECISION HVEC(JDMB) DOUBLE PRECISION DTERM CHARACTER ALGO CHARACTER*30 FNAMEB CHARACTER*80 CBUF CC-----LIST OF VARIABLES WHICH ENTER "SUBLMS" (WITHOUT CHANGING INSIDE) CC NVAR, NCAS, JCST, NDXX, NDXY, NMXV, PREC, CC X, Y, XMED, XMAD, JDIAG, ALGO, FNAMEB CC-----LIST OF VARIABLES THE VALUE OF WHICH IS CALCULATED IN "SUBLMS" CC JQU, JREP, DRCTA, JNDA, RESDU, WEIGHTS, UHYB, AVW, NZWE, RSQ, CC JHYB, JDMB, JERD, JLP, NERR, DTERM CC-----LIST OF VARIABLES WHICH ARE NEEDED FOR THE CALCULATIONS CC C, HVEC, JNDVC, AW CC-----LIST OF SUBROUTINES USED IN "SUBLMS" CC SUBROUTINE SUBREP CC SUBROUTINE GENPN CC SUBROUTINE RANPN CC SUBROUTINE EQUAT CC SUBROUTINE RESTT CC SUBROUTINE RANGS CC SUBROUTINE SHHLF CC FUNCTION PULL CC-----INITIALIZATIONS JERD=0 JFRST=0 JEY=0 AVW=0.0 NZWE=0 JREADZ=0 JHYB=JDIAG AL=NCAS ANVAR=NVAR CC-----CALCULATION OF JQU, THE RANK OF THE ABSOLUTE RESIDUAL TO BE MINIMIZED JQU = INT(AL/2.0) + INT( (ANVAR+1.0)/2.0 ) CC-----CALCULATE THE NUMBER OF REPLICATIONS CALL SUBREP(NVAR,NCAS,NDXX,ALGO,JNDVC,JREP) CC-----MAXIMAL NUMBER OF REPLICATIONS MAXRP=3*JREP NVSB=NVAR-1 NVAD=NVAR+1 CC-----START OF THE ALGORITHM : THE "10" LOOP PERFORMS THE REPLICATIONS DO 10 JJJ=1,JREP CC-----IF (JNDVC(NVAD) = 11 THEN TAKE A COMBINATION OF NVAR POINTS CC-----OUT OF NCAS, ELSE TAKE A RANDOM SAMPLE OF NVAR POINTS OUT OF NCAS 20 IF(JNDVC(NVAD).EQ.11) THEN CALL GENPN(NCAS,NVAR,JNDVC,NDXX,JJJ) ELSE IF (JJJ.GT.2) THEN CALL RANPN(NCAS,NVAR,JNDVC,NDXX,JEY,JLP,MAXRP) IF (JLP.GT.MAXRP) GOTO 210 ELSE IF (JJJ.EQ.2) THEN JLP=JLP+1 DO 30 MNDX=1,NVAR 30 JNDVC(MNDX)=NCAS-MNDX+1 ELSE JLP=JLP+1 DO 40 MNDX=1,NVAR 40 JNDVC(MNDX)=MNDX ENDIF ENDIF ENDIF CC-----STORE THE NVAR SELECTED OBSERVATIONS IN THE MATRIX C DO 50 MNDX=1,NVAR DO 60 KJ=1,NVAD NDR=JNDVC(MNDX) C(MNDX,KJ)=X(KJ,NDR) 60 CONTINUE 50 CONTINUE CC-----THE FOLLOWING FIVE LINES ARE INSTRUCTIONS FOR INFORMING CC-----THE USER ABOUT HOW FAR THE CALCULATIONS OF THE LMS ARE ALREADY ADVANCED CC-----(THESE LINES MAY BE DELETED IF USED ON OTHER COMPUTER) JREADY=(JJJ*100)/JREP IF (MOD(JREADY,5).EQ.0.AND.JREADY.NE.JREADZ) THEN WRITE(CBUF,65) JREADY 65 FORMAT(I4,' PERCENT OF THE CALCULATIONS HAVE BEEN EXECUTED.') CALL STTPUT(CBUF,ISTAT) ENDIF IF (FNAMEB.EQ.'CON'.AND.JREADY.GE.100) PAUSE ' ' JREADZ=JREADY CC-----SOLVE THE SYSTEM OF EQUATIONS STORED IN C IF(NVAR.GT.1) THEN CALL EQUAT(C,NMXV,NDXX,HVEC,JDMB,NVAR,1,NERR,DTERM) IF(NERR.GE.0) GOTO 70 JERD=JERD+1 IF (JNDVC(NVAD).NE.11.AND.JJJ.GT.2) GOTO 20 GOTO 10 ELSE IF(C(1,1).NE.0.0) C(1,1)=C(1,2)/C(1,1) ENDIF CC-----CALCULATE THE RESIDUALS AND THE ASSOCIATED OBJECTIVE FUNCTION CC-----VALUE FOR THE TRIAL ESTIMATE OBTAINED ABOVE 70 DO 80 JNC=1,NCAS RESDU(JNC)=0.0 DO 90 J=1,NVAR RESDU(JNC)=RESDU(JNC)+C(J,1)*X(J,JNC) 90 CONTINUE RESDU(JNC)=ABS(X(NVAD,JNC)-RESDU(JNC)) 80 CONTINUE FCTVA=PULL(AW,NDXY,RESDU,NCAS,JQU)*1.4826 CC-----INITIALIZE THE "BEST" ESTIMATE AND OBJECTIVE FUNCTION VALUE IF(JFRST.NE.1) THEN JFRST=1 DRCTA(NVAD)=FCTVA DO 100 J=1,NVAR JNDA(J)=JNDVC(J) DRCTA(J)=C(J,1) 100 CONTINUE IF (FCTVA.LT.PREC) JHYB=0 CC-----IF DIAGNOSTICS ARE WANTED, THEN CALCULATE THE RDi IF (JHYB.EQ.0) GOTO 200 DO 110 JNC=1,NCAS UHYB(JNC)=RESDU(JNC)/FCTVA*1.4826 110 CONTINUE ELSE 120 IF(FCTVA.LT.PREC) JHYB=0 IF (JHYB.NE.0) THEN DO 130 JNC=1,NCAS UH=RESDU(JNC)/FCTVA*1.4826 IF(UHYB(JNC).LT.UH) UHYB(JNC)=UH 130 CONTINUE ENDIF CC-----UPDATE THE "BEST" ESTIMATE AND OBJECTIVE FUNCTION VALUE CC-----IF NECESSARY IF (FCTVA.LT.DRCTA(NVAD)) THEN DRCTA(NVAD)=FCTVA DO 140 J=1,NVAR JNDA(J)=JNDVC(J) 140 DRCTA(J)=C(J,1) ENDIF ENDIF CC-----IF FCTVA IS LESS THAN PREC (A PRECISION CONSTANT) THEN ONE IS CC-----RESORTED TO THE EXACT FIT CASE 200 CONTINUE C write(6,*) 'FCTVA,PREC', FCTVA, PREC IF (FCTVA.LE.PREC) THEN DRCTA(NVAD)=FCTVA CALL RESTT(DRCTA,1,0,NVAR,NCAS,NVAD,NZWE,X,Y,RESDU,WEIGHTS, 1 XMED,XMAD,NDXX,NDXY,0,SCAL,QUAN,PREC) AVW=(NZWE*1.0 )/AL C RETURN ENDIF 10 CONTINUE 210 CONTINUE CC-----ADJUSTING THE INTERCEPT IN THE CASE OF A MODEL WITH CC-----CONSTANT TERM write(CBUF,*) 'JCST = ',JCST CALL STTPUT(CBUF,ISTAT) IF (JCST.NE.0) THEN CALL RESTT(DRCTA,0,0,NVSB,NCAS,NVAD,NZWE,X,Y,RESDU, 1 WEIGHTS,XMED,XMAD,NDXX,NDXY,0,SCAL,QUAN,PREC) CALL RANGS(RESDU,NCAS) CALL SHHLF(RESDU,NCAS,JQU,SLUTN,BSTD,PREC) DRCTA(NVAR)=SLUTN DRCTA(NVAD)=BSTD*1.4826 ENDIF CC-----COMPUTING THE FINAL SCALE ESTIMATE OF THE LMS AND THE CC-----RESULTING WEIGHTS QUAN=DRCTA(NVAD) QUAN=QUAN*(5.0/(AL-ANVAR)+1.0) CALL RESTT(DRCTA,1,0,NVAR,NCAS,NVAD,NZWE, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NDXX,NDXY,1,SCAL,QUAN,PREC) DRCTA(NVAD)=SCAL AVW=(NZWE*1.0)/AL CC-----COEFFICIENT OF DETERMINATION CORRESPONDING TO THE LMS RSQ= DRCTA(NVAD)*XMAD(NVAD) DO 220 JNC=1,NCAS 220 AW(JNC)=ABS(Y(JNC)*XMAD(NVAD)) RSY=PULL(AW,NDXY,AW,NCAS,JQU)*1.4826 RSQ=1.0-(RSQ/RSY)**2 IF (RSQ.LT.0.0) RSQ=0.0 IF (RSQ.GT.1.0) RSQ=1.0 RETURN END CC ------------------------------------------------------------------------ CC PRINTING LMS ESTIMATES AND SOME OPTIONS CC ------------------------------------------------------------------------ SUBROUTINE PRTLMS(NVAR,NCAS,JCST,NVSB,NVAD,NDXX,NDXY,NMXV,JNDVC, 1 XMED,XMAD,DRCTA,AW,RESDU,WEIGHTS,X,Y,NMVAL,UHYB,AVW,NZWE,RSQ, 1 JREG,JPRT,JPLT,JREP,JERD,JLP,JQU,JHYB,JDIAG,NSTOP,JHEAD,LAB, 1 PREC,LUB) INCLUDE 'MID_REL_INCL:implicit.inc' DIMENSION X(NDXX,NDXY),Y(NDXY),XMED(NDXX),XMAD(NDXX),UHYB(NDXY) DIMENSION DRCTA(NMXV),AW(NDXY),RESDU(NDXY),NMVAL(NDXY) DIMENSION WEIGHTS(NDXY) DIMENSION JNDVC(NDXX) CHARACTER*10 LAB(NDXX) CHARACTER*60 JHEAD WRITE(LUB,8000) WRITE(LUB,8010) IF (JPRT.NE.0) WRITE(LUB,8020) JQU NSTOP=0 IF (JERD.NE.0) THEN IF (JNDVC(NVAD).NE.11) THEN IF (JERD.EQ.JLP) THEN WRITE(LUB,8030) NVAR,NCAS NSTOP=1 RETURN ELSE IF (JPRT.NE.0) WRITE(LUB,8040) JLP,NVAR,NCAS,JERD ENDIF ELSE IF (JERD.EQ.JREP) THEN WRITE(LUB,8050) NVAR,NCAS NSTOP=1 RETURN ELSE IF (JPRT.NE.0) WRITE(LUB,8060) JERD,NVAR,NCAS,JREP ENDIF ENDIF ENDIF CC-----LMS IN THE CASE OF EXACT FIT IF (DRCTA(NVAD).LE.PREC) THEN WRITE(LUB,8070) DO 10 J=1,NVSB 10 WRITE(LUB,8080) LAB(J),DRCTA(J) IF (JCST.EQ.0) WRITE(LUB,8080) LAB(NVAR),DRCTA(NVAR) IF (JCST.NE.0) WRITE(LUB,8090) DRCTA(NVAR) WRITE(LUB,8100) DRCTA(NVAD) WRITE(LUB,8110) NZWE JREG=2 CC-----RESIDUAL PLOT FOR THE LMS IN THE CASE OF EXACT FIT CALL RDUAL(DRCTA,0,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC, 1 JREG, X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL, 1 NDXX,NDXY,JHEAD,LAB) NSTOP=1 RETURN ENDIF CC-----LMS IN THE GENERAL CASE WRITE(LUB,8130) DO 20 J=1,NVSB 20 WRITE(LUB,8080) LAB(J),DRCTA(J) IF (JCST.EQ.0) WRITE(LUB,8080) LAB(NVAR),DRCTA(NVAR) IF (JCST.NE.0) WRITE(LUB,8090) DRCTA(NVAR) WRITE(LUB,8200) DRCTA(NVAD) CC-----PRINT COEFFICIENT OF DETERMINATION OF THE LMS WRITE(LUB,8210) RSQ JREG=2 CC-----RESIDUAL PLOTS IF REQUESTED CALL RDUAL(DRCTA,1,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL,NDXX,NDXY, 1 JHEAD,LAB) IF (JHYB.NE.0.AND.JDIAG.EQ.1) THEN WRITE(LUB,8220) UH=AMDAN(AW,NDXY,UHYB,NCAS) DO 40 JNC=1,NCAS RDI=UHYB(JNC)/UH WRITE(LUB,8230) NMVAL(JNC),UHYB(JNC),RDI 40 CONTINUE ENDIF 8000 FORMAT(/1X,78('*')/) 8010 FORMAT(' LEAST MEDIAN OF SQUARES REGRESSION'/ 1 1X,34('*')/) 8020 FORMAT(' THE MINIMIZATION OF THE ',I4, 1 'TH ORDERED SQUARED RESIDUAL IS PERFORMED.') 8030 FORMAT(//' ALL SUBSAMPLES OF ',I3,' POINTS OUT OF ',I5, 1 ' LED TO A SINGULAR SYSTEM OF EQUATIONS.'/ 1 ' THIS MEANS THAT THERE IS AN ENORMOUS PROBLEM OF', 1 ' COLLINEARITY IN THE DATA SET.'/ 1 ' TRY THE ANALYSIS AGAIN IN A LOWER DIMENSION.'//) 8040 FORMAT(/' ON A TOTAL OF ',I8,' SUBSAMPLES (OF', 1 I3,'POINTS OUT OF ',I5,')', /I7, 1 ' SUBSAMPLES LED TO A SINGULAR SYSTEM OF EQUATIONS.'/ 1 ' THE SOLUTION IS ONLY BASED ON THE GOOD SUBSAMPLES.'//) 8050 FORMAT(//' ALL COMBINATIONS OF ,I3,15H POINTS OUT OF ',I5, 1 ' LED TO A SINGULAR SYSTEM OF EQUATIONS.'/ 1 ' THIS MEANS THAT THERE IS AN ENORMOUS PROBLEM OF', 1 ' COLLINEARITY IN THE DATA SET.'/ 1 ' TRY THE ANALYSIS AGAIN IN A LOWER DIMENSION.'//) 8060 FORMAT(/' THERE WERE ',I7,' COMBINATIONS (OF ',I3, 1 ' POINTS OUT OF ,I5,16H) OF A TOTAL OF '/I8,' COMBINATIONS,', 1 ' WHICH LED TO A SINGULAR SYSTEM OF EQUATIONS.'/ 1 ' THE SOLUTION IS ONLY BASED ON THE GOOD COMBINATIONS.'//) 8070 FORMAT(/' MORE THAN THE HALF OF THE DATA LIE ON THE' 1 /' SAME HYPERPLANE, WITH COEFFICIENTS = '//) 8080 FORMAT(3X,A10,F20.5) 8090 FORMAT(3X,' CONSTANT',F20.5) 8100 FORMAT(/' SCALE ESTIMATE',10X,'= ',F15.5/) 8110 FORMAT(' THERE ARE',2X,I4,' POINTS ON THE HYPERPLANE') 8130 FORMAT(/5X,'VARIABLE',9X,'COEFFICIENT'/3X,30('-')) 8200 FORMAT(/' FINAL SCALE ESTIMATE',9X,'= ',F15.5) 8210 FORMAT(/' COEFFICIENT OF DETERMINATION = ',F15.5/) 8220 FORMAT(//1X,'INDEX',3X,' Ui ',3X,'RESISTANT '/ 1 25X,'DIAGNOSTIC'/1X,34('-')) 8230 FORMAT(1X,I4,3X,F13.4,3X,F10.4) RETURN END