CC$large CC CC PROGRESS : PROGRAM FOR ROBUST REGRESSION CC----------------------------------------------------------------------------- CC-----LIST OF VARIABLES AND ARRAYS IN THIS PROGRAM CC----------------------------------------------------------------------------- CC MAXN = maximal number of cases CC MAXP = maximal number of explanatory variables CC MAXP1 = MAXP+1 = maximal number of coefficients, CC including intercept CC MAXPP1 = MAXP*(MAXP+1) CC----------------------------------------------------------------------------- CC X(MAXP1,MAXN) = matrix of the values of the explanatory CC variables (transposed design matrix) CC Y(MAXN) = vector with the values of the response variable CC RESDU(MAXN) = vector for storing the residuals CC WEIGHTS(MAXN) = weight assigned to each observation CC AA(MAXN) = working vector CC AM(MAXN) = working vector CC AW(MAXN) = working vector CC AW2(MAXN) = working vector CC A(MAXP+1) = vector of regression coefficients (ls or rls) CC DA(MAXP+1) = vector with the stand. deviation of the coefficients CC COEFF(MAXP+1) = vector of coefficients for the lms/lts CC XMED(MAXP+1) = vector with the medians of all variables CC XMAD(MAXP+1) = vector with median absolute deviations CC UHYB(MAXN) = vector with the resistant diagnostics CC C(MAXP,MAXP+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(MAXP,MAXP+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(MAXPP1) = auxiliary vector CC----------------------------------------------------------------------------- CC JNDEX(MAXN) = vector with the original index of the cases CC JSUBS(MAXP+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(MAXP+1) = vector with the "missing value codes" CC JMISS(MAXP+1) = JMISS(J), (J=1,NVAR+1) is 1 if var(j) contains missing CC values and is 0 otherwise CC JBEST(MAXP) = vector containing the index of the cases CC which yield the "best" lms/lts solution CC JPLACE(MAXP+1) = vector containing the positions of the variables CC used in the analysis (from the total data set) CC LAB(MAXP+1) = vector with the label of each variable CC JREPTAB(MAXP+1)= default number of random subsets (depending on CC number of regressors) CC JREPLOW(MAXP+1)= number of observations which allow to take all CC combinations (depending on number of regressors) CC JREPHI(MAXP+1) = number of observations which do not allow to CC take all combinations CC FACLMS(11) = asympt. consistency factors for lms scale estimation CC FACLTS(11) = asympt. consistency factors for lts scale estimation 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 NFAC = NVAR-1 CC NVAD = NVAR+1 CC NCAS = number of cases CC AL = real value of ncas CC ALGO = character which is 'A' for drawing all subsets CC and 'R' for drawing random subsets 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 PREC = precision of the calculations CC JERD = number of subsets giving singular system of equations CC JEY,JLP = integer variables used for counting the number of CC replications in the lms/lts algorithm CC JHYB = 0 if half of the observations lie on a hyperplane CC JPLXY = option for xy plot (0=no xy plot, 1=otherwise) CC NREP = number of (random) subsets to be considered CC JREG = which regression ? (1=ls, 2=lms/lts, 3=rls) CC MROB = which robust method ? (1=lms, 2=lts) CC JDIAG = option for diagnostics (0=no diagnostics, else 1) CC JQU = the quantile to be minimized CC JEDAUL = default quantile which yields highest breakdown value CC NSTOP = 1 means stop the analysis, otherwise NSTOP=0 CC LOCA = 1 for univariate location/scale estimation CC BSTD = univariate lms/lts scale estimate CC SLUTN = univariate lms/lts location estimate CC FACTOR = asympt. consistency factor for lms/lts scale CC FINITF = finite sample correction factor for lms/lts scale CC RDI = robust distance 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 PRTRLS :prints the rls (or ls) results CC SUBROUTINE SUBLMS :algorithm for calculating the lms/lts estimates CC SUBROUTINE PRTLMS :prints the lms/lts 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/lts 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/lts 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 LMSLOC :calculates univariate lms CC SUBROUTINE LTSLOC :calculates univariate lts 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 NCOMB :computes number of combinations of k 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 univariate location and scale estimators CC SUBROUTINE SMISLOC :treatment of missing values in the univariate case 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 SUBROUTINE JQUANT :asks for quantile to be minimized CC FUNCTION NBREAK :computes breakdown value of lms/lts CC----------------------------------------------------------------------------- CC 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 LINES: CC ------------------------------------------------------------------------ DIMENSION X(15,1000),Y(1000),RESDU(1000),UHYB(1000) DIMENSION AW(1000),AW2(1000),WEIGHTS(1000) DIMENSION COEFF(15),XMED(15),XMAD(15),C(14,15) DIMENSION VALMS(15),A(15),DA(15) DOUBLE PRECISION H(14,15),HVEC(210) DIMENSION FACLMS(11),FACLTS(11) INTEGER JNDEX(1000) INTEGER JSUBS(15) INTEGER JBEST(14) INTEGER JPLACE(15) INTEGER JMISS(15) CHARACTER*10 LAB(15) INTEGER JREPTAB(15) INTEGER JREPLOW(15) INTEGER JREPHI(15) CC-----THE FOLLOWING THREE VARIABLES NEED TO HAVE LENGTH = MAXP+1 CC-----YOU DO NOT HAVE TO CHANGE FACLMS AND FACLTS DATA JREPTAB/500,1000,1500,2000,2500,10*3000/ DATA JREPLOW/500,50,22,17,15,14,9*0/ DATA JREPHI/1000000,1414,182,71,43,32,27,24,23,22,22,22,22,23,23/ DATA FACLMS/1.4826,1.3998,1.3238,1.2535,1.1882,1.1272,1.0700, 1 1.0160,0.9648,0.9161,0.8693/ DATA FACLTS/2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660, 1 1.879,1.7973,1.7203,1.6473/ MAXP=14 MAXN=1000 CC ------------------------------------------------------------------- CC INITIALIZATION CC ------------------------------------------------------------------------ LUA=1 LUB=2 LUC=3 MAXP1=MAXP+1 MAXPP1=MAXP*MAXP1 NZWE=0 PREC=1.0E-6 JLP=0 JREG=0 NSTOP=0 DO 5 J=1,MAXP JBEST(J)=0 5 CONTINUE JHEAD=' ' FNAMEA=' ' FNAMEB=' ' FNAMEC=' ' CC ------------------------------------------------------------------------ CC READ THE DATA CC ------------------------------------------------------------------------ CALL RDAT(NCAS,NVAR,JCST,JPRT,NFAC,NVAD,VALMS, 1 X,Y,JPLT,JMISS,PREC,JPLACE,LAB,XMED,XMAD,AW,AW2,JNDEX, 1 RESDU,WEIGHTS,MAXP1,MAXN,MVAL,MAXP,LUA,LUB,LUC,JHEAD, 1 YNSAVE,JDIAG,NSTOP,FNAMEA,FNAMEB,FNAMEC,MROB,FACLMS,FACLTS) IF (NSTOP.EQ.1) GOTO 30 CC ------------------------------------------------------------------------ CC CREATE THE BEGINNING OF THE OUTPUT: THE SELECTED OPTIONS CC ------------------------------------------------------------------------ CALL WRDATA(NVAR,NCAS,JCST,MVAL,JPRT,YNSAVE,FNAMEA,FNAMEB, 1 FNAMEC,JHEAD,LUB,MROB,LAB,MAXP1,NVAD,JPLT) CC ------------------------------------------------------------------------ CC IF THE DATA CONTAIN MISSING VALUES, THEN TREAT THEM CC ------------------------------------------------------------------------ IF (MVAL.NE.0) THEN CALL SMISSING(NVAR,NCAS,MAXP1,MAXN,JCST,X,Y,AW,JNDEX,MVAL, 1 JMISS,JSUBS,VALMS,NSTOP,LAB,JPRT,LUB) IF (NSTOP.EQ.1) GOTO 20 ENDIF CC ------------------------------------------------------------------------ CC ASK FOR QUANTILE AND NUMBER OF SUBSETS USED IN ROBUST REGRESSION CC ------------------------------------------------------------------------ CALL JQUANT(NCAS,NVAR,MROB,JQU,FACTOR,FACLMS,FACLTS) CALL SUBREP(NCAS,NVAR,MAXP1,ALGO,NREP,JREPTAB,JREPLOW,JREPHI) CC ------------------------------------------------------------------------ CC CREATE SECOND PART OF THE OUTPUT: THE DATA AND AN X-Y PLOT CC ------------------------------------------------------------------------ CALL WRDATB(NVAR,NCAS,MAXP1,MAXN,JCST,MVAL,LAB,JNDEX,X,Y, 1 AW,JPLT,JPRT,JREG,JHEAD,LUB,MROB) CC ------------------------------------------------------------------------ CC STANDARDIZATION OF THE DATA (WITH OPTIONAL PRINTING) CC AND CALCULATION OF SOME STATISTICS ON THE VARIABLES CC ------------------------------------------------------------------------ CALL STATIS(X,Y,XMED,XMAD,AW2,WEIGHTS,AW,JNDEX,C,H,JSUBS, 1 NSTOP,NVAR,NVAD,NFAC,NCAS,JCST,JPRT,MAXP1,MAXN,LUB,MAXP, 1 PREC,LAB) IF (NSTOP.EQ.1) GOTO 20 CC ------------------------------------------------------------------------ CC CALCULATION OF THE LS ESTIMATES CC ------------------------------------------------------------------------ IF (NCAS.LE.(NVAR*2)) THEN WRITE(*,8050) NCAS,NVAR IF (JCST.EQ.1) WRITE(*,8060) GOTO 20 ENDIF JREG=1 IF (NVAR.EQ.1) THEN CALL LSL(NCAS,MAXN,X,Y,WEIGHTS,A,FCKW,H,MAXP,MAXP1) ELSE CALL LSREG(MAXP1,MAXN,MAXP,NVAR,NCAS,NVAR,A,X,Y,WEIGHTS,DA, 1 H,FCKW,HVEC,MAXPP1,JMISS) ENDIF CC ------------------------------------------------------------------------ CC RETRANSFORMATION OF THE LS ESTIMATES CC ------------------------------------------------------------------------ CALL RTRAN(NVAR,JCST,NFAC,NVAD,MAXP1,XMED,XMAD,A,NVAD,FCKW) CALL TRC(H,DA,MAXP,MAXP1,NVAR,JCST,NFAC,NVAD,XMED,XMAD) CC ------------------------------------------------------------------------ CC PRINTING THE LS RESULTS WITH THE REQUESTED OPTIONS CC ------------------------------------------------------------------------ CALL PRTRLS(NVAR,NCAS,JCST,NFAC,NVAD,MAXP1,MAXN,MAXP, 1 XMED,XMAD,A,DA,H,AW,RESDU,WEIGHTS,X,Y,JNDEX,AVW,NZWE,FCKW, 1 JREG,JPRT,JPLT,JHEAD,LAB,PREC,LUB,MROB) CC ------------------------------------------------------------------------ CC CALCULATION OF LMS REGRESSION CC ------------------------------------------------------------------------ WRITE(*,8000) CALL SUBLMS(NVAR,NCAS,JCST,MAXP1,MAXN,MAXP,PREC,X,Y,XMED,XMAD, 1 JDIAG,ALGO,FNAMEB,JQU,NREP,COEFF,JBEST,RESDU,WEIGHTS,UHYB,AVW, 1 NZWE,RSQ,JHYB,MAXPP1,JERD,JLP,C,HVEC,JSUBS,AW,AW2, 1 MROB,FCKW,PRELSCAL,FACTOR) CC ------------------------------------------------------------------------ CC RETRANSFORMATION OF THE LMS ESTIMATES CC ------------------------------------------------------------------------ CALL RTRAN(NVAR,JCST,NFAC,NVAD,MAXP1,XMED,XMAD,COEFF,NVAD,FCKW) COEFF(NVAD)=COEFF(NVAD)*XMAD(NVAD) FCKW=FCKW/XMAD(NVAD) PRELSCAL=PRELSCAL*XMAD(NVAD) CC ------------------------------------------------------------------------ CC PRINTING THE LMS RESULTS WITH THE REQUESTED OPTIONS CC ------------------------------------------------------------------------ CALL PRTLMS(NVAR,NCAS,JCST,NFAC,NVAD,MAXP1,MAXN,MAXP,XMED, 1 XMAD,COEFF,AW,RESDU,WEIGHTS,X,Y,JNDEX,UHYB,NZWE,RSQ,JREG,JBEST, 1 JPRT,JPLT,NREP,JERD,JLP,JQU,JHYB,JDIAG,NSTOP,JHEAD,LAB,PREC,LUB, 1 FCKW,MROB,ALGO,PRELSCAL) IF (NSTOP.NE.1) THEN CC ------------------------------------------------------------------------ CC CALCULATION OF RLS REGRESSION CC ------------------------------------------------------------------------ JREG=3 WRITE(*,8005) DO 10 J=1,MAXP1 A(J)=0.0 DA(J)=0.0 DO 10 JNC=1,MAXP 10 H(JNC,J)=0.D0 IF (NVAR.EQ.1) THEN CALL LSL(NCAS,MAXN,X,Y,WEIGHTS,A,FCKW,H,MAXP,MAXP1) ELSE CALL LSREG(MAXP1,MAXN,MAXP,NVAR,NCAS,NVAR,A,X,Y,WEIGHTS,DA, 1 H,FCKW,HVEC,MAXPP1,JMISS) ENDIF CC ------------------------------------------------------------------------ CC RETRANSFORMATION OF THE RLS ESTIMATES CC ------------------------------------------------------------------------ CALL RTRAN(NVAR,JCST,NFAC,NVAD,MAXP1,XMED,XMAD,A,NVAD,FCKW) CALL TRC(H,DA,MAXP,MAXP1,NVAR,JCST,NFAC,NVAD,XMED,XMAD) CC ------------------------------------------------------------------------ CC PRINTING THE RLS RESULTS WITH THE REQUESTED OPTIONS CC ------------------------------------------------------------------------ CALL PRTRLS(NVAR,NCAS,JCST,NFAC,NVAD,MAXP1,MAXN,MAXP, 1 XMED,XMAD,A,DA,H,AW,RESDU,WEIGHTS,X,Y,JNDEX,AVW,NZWE,FCKW, 1 JREG,JPRT,JPLT,JHEAD,LAB,PREC,LUB,MROB) ENDIF CC END OF THE RUN WRITE (*,8010) 20 IF (YNSAVE.EQ.'Y') WRITE (*,8020) FNAMEC IF (FNAMEA.NE.'CON') WRITE(*,8030) FNAMEA IF (FNAMEB.NE.'CON'.AND.FNAMEB.NE.'PRN') WRITE(*,8040)FNAMEB CC ----------------------------------------------------------------- CC FORMATS CC ----------------------------------------------------------------- 8000 FORMAT(//' Algorithm for calculating the positive breakdown ', 1 'regression is started.'//) 8005 FORMAT(//45h Progress is presently calculating reweighted, 1 18h least squares .../) 8010 FORMAT(//41h The run has been executed successfully. //) 8020 FORMAT(' The data are saved in the file: ',A30) 8030 FORMAT(' The data has been read from the file: ',A30) 8040 FORMAT(' The output has been written in the file: ',A30) 8050 FORMAT(//46h Too many coefficients according to the number, 1 10h of cases.//16h Number of cases,8X,2H= ,I5/ 1 26h Number of coefficients = ,I5,/25h The number of cases must, 1 37h be twice the number of coefficients.) 8060 FORMAT(31h (Including the constant term!)) 30 STOP 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,YNSAVE,FNAMEA, 1 FNAMEB,FNAMEC,JHEAD,LUB,MROB,LAB,MAXP1,NVAD,JPLT) CHARACTER YNSAVE CHARACTER*30 FNAMEA,FNAMEB,FNAMEC CHARACTER*60 JHEAD CHARACTER*10 LAB(MAXP1) WRITE(LUB,8000) WRITE(LUB,8005) WRITE(LUB,8006) WRITE(LUB,8010) JHEAD IF (JCST.EQ.0) THEN WRITE(LUB,8020) NCAS,NVAR,LAB(NVAD) WRITE(LUB,8040) ELSE WRITE(LUB,8020) NCAS,NVAR-1,LAB(NVAD) WRITE(LUB,8030) ENDIF IF (MROB.EQ.1) WRITE(LUB,8050) IF (MROB.EQ.2) WRITE(LUB,8060) IF (FNAMEB.EQ.'CON') PAUSE ' ' IF (YNSAVE.EQ.'Y') WRITE(LUB,8080) FNAMEC IF (FNAMEA.NE.'CON') WRITE(LUB,8070) FNAMEA WRITE(LUB,8085) FNAMEB IF (JPRT.EQ.0) WRITE(LUB,8090) IF (JPRT.EQ.1) WRITE(LUB,8100) IF (JPRT.EQ.2) WRITE(LUB,8110) IF (JPLT.EQ.0) WRITE(LUB,8200) IF (JPLT.EQ.1) WRITE(LUB,8210) IF (JPLT.EQ.2) WRITE(LUB,8220) IF (JPLT.EQ.3) WRITE(LUB,8230) IF (MVAL.EQ.0) WRITE(LUB,8120) IF (MVAL.EQ.1) WRITE(LUB,8130) IF (MVAL.EQ.2) WRITE(LUB,8140) 8000 FORMAT(//30X,19(1H*)/30X,19H* P R O G R E S S */ 1 30X,19(1H*)///, 1 ' This program performs a robust regression analysis using', 1 ' either'/ 1 ' the least median of squares (LMS) method '/ 1 ' or'/ 1 ' the least trimmed squares (LTS) method'/ 1 ' both described in :'// 1 ' P.J. 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.J. Rousseeuw and A.M. Leroy (1987), Robust Regression'/ 1 ' and Outlier Detection, Wiley, New York.'//) 8005 FORMAT(' Since January 16, 1996 the program was upgraded'/ 1 ' in collaboration with Mia Hubert, yielding:'/ 1 ' - better optimization by adjusting the intercept at each', 1 ' slope '/ 1 ' (this yields the exact LMS for simple regression)') 8006 FORMAT(' - optimal default value: h = [(n+p+1)/2]'/ 1 ' - option to set a lower breakdown value (higher h) '/ 1 ' - choice between exhaustive and random p-subsets '/ 1 ' - choice of the number of random p-subsets '/ 1 ' - new definition of robust R-squared.'/ 1 ' Therefore, the outputs differ slightly from those', 1 ' in the 1987 book.'/' More information on these recent', 1 ' changes can be found in :'// 1 ' P.J. Rousseeuw and M. Hubert (1997), Recent developments'/ 1 ' in PROGRESS, in L1-Statistical Procedures and Related'/ 1 ' Topics, edited by Y. Dodge. Institute of Mathematical'/ 1 ' Statistics Lecture Notes-Monograph Series, Volume 31,'/ 1 ' Hayward, California, pages 201-214.'//) 8010 FORMAT(/' Title: ',A60) 8020 FORMAT(' Number of cases in the data set:',19X,I5/ 1 ' Number of explanatory variables in the regression: ',I5/ 1 ' The response variable is: ',A10) 8030 FORMAT(' You have chosen a model with a constant term', 1 ' (regression with intercept).') 8040 FORMAT(' You have chosen a model without constant term', 1 ' (no intercept).') 8050 FORMAT(' You have chosen the least median of squares', 1 ' (LMS) method.') 8060 FORMAT(' You have chosen the least trimmed squares', 1 ' (LTS) method.') 8070 FORMAT(' Your data reside in the file: ',A30) 8080 FORMAT(' The data are saved in the file: ',A30) 8085 FORMAT(' This output is written in the file: ',A30) 8090 FORMAT(' You have chosen small output.') 8100 FORMAT(' You have chosen medium-sized output.') 8110 FORMAT(' You have chosen large output.') 8200 FORMAT(' You have chosen not to plot residuals.') 8210 FORMAT(' You have chosen a plot of residuals versus ', 1 'estimated y.') 8220 FORMAT(' You have chosen a plot of residuals versus ', 1 'their index.') 8230 FORMAT(' You have chosen a plot of residuals versus their index', 1 ' and '/ 1 ' versus estimated y.') 8120 FORMAT(' The data is assumed not to have missing values.') 8130 FORMAT(' Treatment of missing values in option 1:'/ 1 ' cases with at least one missing value will be deleted.'/) 8140 FORMAT(' Treatment of missing values in option 2:'/ 1 ' First, any case with a missing value for the response', 1 ' variable'/ 1 ' or for all explanatory variables will be deleted.'/ 1 ' For the remaining cases, each missing value is replaced by'/ 1 ' the median of the non-missing values.'/) 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,MAXP1,MAXN,JCST,MVAL,LAB,JNDEX,X,Y, 1 AW,JPLT,JPRT,JREG,JHEAD,LUB,MROB) DIMENSION X(MAXP1,MAXN),Y(MAXN),AW(MAXN) INTEGER JNDEX(MAXN) CHARACTER*10 LAB(MAXP1) CHARACTER*60 JHEAD NFAC=NVAR-1 NVAD=NVAR+1 JPLXY=0 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,NFAC),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) JNDEX(JNC),(X(J,JNC),J=1,NFAC), 1 X(NVAD,JNC) IF (JCST.EQ.0) WRITE(LUB,8140) JNDEX(JNC),(X(J,JNC),J=1,NVAD) GOTO 5 15 WRITE(LUB,8140) JNDEX(JNC),(X(J,JNC),J=1,6) IF (JCST.EQ.1.AND.NFAC.GE.7) 1 WRITE(LUB,8150) (X(J,JNC),J=7,NFAC),X(NVAD,JNC) IF (JCST.EQ.1.AND.NFAC.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,JNDEX,MAXN,LUB,JREG, 1 JHEAD,MAXP1,NVAD,LAB,MROB) 8110 FORMAT(/22h The observations are:/) 8120 FORMAT(/' The remaining observations are: '/) 8130 FORMAT(' i ',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 PRINTING THE RLS ESTIMATES AND SOME OPTIONS CC ------------------------------------------------------------------------ SUBROUTINE PRTRLS(NVAR,NCAS,JCST,NFAC,NVAD,MAXP1,MAXN,MAXP, 1 XMED,XMAD,A,DA,H,AW,RESDU,WEIGHTS,X,Y,JNDEX,AVW,NZWE,FCKW, 1 JREG,JPRT,JPLT,JHEAD,LAB,PREC,LUB,MROB) DIMENSION X(MAXP1,MAXN),Y(MAXN),WEIGHTS(MAXN) DIMENSION XMED(MAXP1),XMAD(MAXP1),AW(MAXN),RESDU(MAXN) DIMENSION A(MAXP),DA(MAXP) DOUBLE PRECISION H(MAXP,MAXP1),DS,PTVAL,DST,PVAL INTEGER JNDEX(MAXN) CHARACTER*10 LAB(MAXP1) CHARACTER*60 JHEAD AL=NCAS ANVAR=NVAR WRITE(LUB,8000) IF (JREG.EQ.1) THEN WRITE(LUB,8010) NNUL=NCAS AVW=1 ELSE IF (MROB.EQ.1) THEN WRITE(LUB,8020) ELSE WRITE(LUB,8025) ENDIF NNUL=NZWE ENDIF NDF=NNUL-NVAR CC-----PRINTING THE COEFFICIENTS WITH THEIR STANDARD ERROR, CC-----AND T-VALUE WRITE(LUB,8030) DO 10 J=1,NFAC 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 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 CC-----CALCULATION OF THE COEFFICIENT OF DETERMINATION AND F-VALUE RSQ=RSQU(NCAS,NVAD,JCST,Y,MAXN,FCKW,FVALUE,PREC, 1 XMAD,XMED,MAXP1,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 IF (JREG.EQ.1) THEN WRITE(LUB,8070) A(NVAD) ELSE WRITE(LUB,8075) A(NVAD) ENDIF CC-----PRINTING THE VARIANCE-COVARIANCE MATRIX CALL SCHCV(NVAR,MAXP1,MAXP,H,LUB) FCKW=FCKW/(AL*AVW-ANVAR) CC-----PRINT COEFFICIENT OF DETERMINATION, F-VALUE, AND P-VALUE IF (JREG.EQ.1) THEN WRITE(LUB,8100) RSQ ELSE WRITE(LUB,8105) RSQ ENDIF 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,JPLT,AW,JNDEX,MAXP1,MAXN, 1 JHEAD,LAB,MROB) ENDIF 8000 FORMAT(/1X,78(1H*)/) 8010 FORMAT(25H LEAST SQUARES REGRESSION /1X,24(1H*)/) 8020 FORMAT(42H REWEIGHTED LEAST SQUARES BASED ON THE LMS/ 1 1X,41(1H*)/) 8025 FORMAT(42H REWEIGHTED LEAST SQUARES BASED ON THE LTS/ 1 1X,41(1H*)/) 8030 FORMAT(/5X,8hvariable,5X,11hcoefficient, 1 4X,12hstand. error,5X,'t-statistic',5X,'p-value'/ 1 3X,70(1H-)) 8040 FORMAT(3X,A10,5X,F11.5,4X,F12.5,7X,F9.5,3X,F9.5) 8050 FORMAT(5X,8hconstant,5X,F11.5,4X,F12.5,7X,F9.5,3X,F9.5) 8060 FORMAT(/20h Degrees of freedom ,5X,2H= ,4X,I5) 8070 FORMAT(/' LS scale estimate',7X,2H= ,F15.5) 8075 FORMAT(/' RLS scale estimate',6X,2H= ,F15.5) 8080 FORMAT(/15h Sum of squares,10X,2H= ,F15.5) 8090 FORMAT(/24h Weighted sum of squares,3H = ,F15.5) 8100 FORMAT(/' R-squared (coefficient of determination) = ',F7.5/) 8105 FORMAT(/' Weighted R-squared (coefficient of determination) = ', 1 F7.5/) 8110 FORMAT(' F-statistic = ',F12.3,7h (with ,I3,5h and , 1 I4,4h df),5X,'p-value = ',F7.5) 8200 FORMAT(/' There are ',I4,' points with nonzero weight.') 8210 FORMAT(/' Average weight = ',F7.5/) RETURN END CC ------------------------------------------------------------------------ CC SUBROUTINE FOR CALCULATING THE LMS/LTS REGRESSION CC ------------------------------------------------------------------------ SUBROUTINE SUBLMS(NVAR,NCAS,JCST,MAXP1,MAXN,MAXP,PREC,X,Y,XMED, 1 XMAD,JDIAG,ALGO,FNAMEB,JQU,NREP,COEFF,JBEST,RESDU,WEIGHTS,UHYB, 1 AVW,NZWE,RSQ,JHYB,MAXPP1,JERD,JLP,C,HVEC,JSUBS,AW,AW2, 1 MROB,FCKW,PRELSCAL,FACTOR) DIMENSION X(MAXP1,MAXN),Y(MAXN),C(MAXP,MAXP1),COEFF(MAXP1) DIMENSION XMED(MAXP1),XMAD(MAXP1) DIMENSION AW(MAXN),RESDU(MAXN),UHYB(MAXN),WEIGHTS(MAXN),AW2(MAXN) DOUBLE PRECISION HVEC(MAXPP1) INTEGER JBEST(MAXP),JSUBS(MAXP1) CHARACTER ALGO CHARACTER*30 FNAMEB CC-----LIST OF VARIABLES WHICH ENTER "SUBLMS" (WITHOUT CHANGING INSIDE) CC NVAR, NCAS, JCST, MAXP1, MAXN, MAXP, MAXPP1, PREC, CC X, Y, XMED, XMAD, JDIAG, ALGO, FNAMEB, JQU, NREP, MROB, FACTOR CC-----LIST OF VARIABLES THE VALUE OF WHICH IS CALCULATED IN "SUBLMS" CC COEFF, JBEST, RESDU, WEIGHTS, UHYB, AVW, NZWE, RSQ, CC JHYB, JERD, JLP, FCKW, PRELSCAL CC-----LIST OF VARIABLES WHICH ARE NEEDED FOR THE CALCULATIONS CC C, HVEC, JSUBS, AW, AW2 CC-----LIST OF SUBROUTINES USED IN "SUBLMS" CC SUBROUTINE GENPN CC SUBROUTINE RANPN CC SUBROUTINE EQUAT CC SUBROUTINE RESTT CC SUBROUTINE RANGS CC SUBROUTINE LMSLOC CC SUBROUTINE LTSLOC CC FUNCTION PULL CC-----INITIALIZATIONS JERD=0 JFRST=0 JEY=0 AVW=0.0 NZWE=0 JREADZ=0 JHYB=1 AL=NCAS ANVAR=NVAR CC-----MAXIMAL NUMBER OF REPLICATIONS MAXRP=3*NREP NFAC=NVAR-1 NVAD=NVAR+1 CC-----START OF THE ALGORITHM : THE "10" LOOP PERFORMS THE REPLICATIONS DO 10 JJJ=1,NREP 20 IF (ALGO.EQ.'A') THEN IF (JJJ.EQ.1) THEN DO 5 J=1,NFAC JSUBS(J)=J 5 CONTINUE JSUBS(NVAR)=NVAR-1 ENDIF CALL GENPN(NCAS,NVAR,JSUBS) JLP=JLP+1 ELSE IF (JJJ.GT.2) THEN IF (JLP.GE.MAXRP) GOTO 210 CALL RANPN(NCAS,NVAR,JSUBS,MAXP1,JEY,JLP,MAXRP) ELSE IF (JJJ.EQ.2) THEN JLP=JLP+1 DO 30 MNDX=1,NVAR 30 JSUBS(MNDX)=NCAS-MNDX+1 ELSE JLP=JLP+1 DO 40 MNDX=1,NVAR 40 JSUBS(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=JSUBS(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 ARE ALREADY ADVANCED CC-----(THESE LINES MAY BE DELETED IF USED ON OTHER COMPUTER) JREADY=(JJJ*100)/NREP IF (MOD(JREADY,5).EQ.0.AND.JREADY.NE.JREADZ) 1 WRITE(*,8010) JREADY 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,MAXP,MAXP1,HVEC,MAXPP1,NVAR,1,NERR) IF(NERR.GE.0) GOTO 71 JERD=JERD+1 IF (ALGO.EQ.'R'.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 71 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)=X(NVAD,JNC)-RESDU(JNC) AW(JNC)=RESDU(JNC) 80 CONTINUE CC-----ADJUST INTERCEPT AND RECALCULATE THE RESIDUALS IF (JCST.EQ.1) THEN CALL RANGS(AW,NCAS) IF (MROB.EQ.1) THEN CALL LMSLOC(AW,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) ELSE CALL LTSLOC(AW,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) ENDIF FCTW=BSTD/FACTOR C(NVAR,1)=C(NVAR,1)+SLUTN DO 85 JNC=1,NCAS RESDU(JNC)=RESDU(JNC)-SLUTN 85 CONTINUE GOTO 99 ENDIF CC-----CALCULATE THE ASSOCIATED OBJECTIVE FUNCTION VALUE CC-----FOR THE TRIAL ESTIMATE OBTAINED ABOVE IF (MROB.EQ.1) THEN DO 95 JNC=1,NCAS AW2(JNC)=ABS(RESDU(JNC)) 95 CONTINUE FCTW=PULL(AW,MAXN,AW2,NCAS,JQU) ELSE DO 97 JNC=1,NCAS AW2(JNC)=RESDU(JNC)**2 97 CONTINUE CALL RANGS(AW2,NCAS) FCTW=0 DO 98 J=1,JQU FCTW=FCTW+AW2(J) 98 CONTINUE FCTW=SQRT(FCTW/JQU) ENDIF 99 IF (FCTW.LT.PREC) THEN JHYB=0 GOTO 200 ENDIF CC-----INITIALIZE THE "BEST" ESTIMATE AND OBJECTIVE FUNCTION VALUE IF(JFRST.NE.1) THEN JFRST=1 COEFF(NVAD)=FCTW DO 100 J=1,NVAR JBEST(J)=JSUBS(J) COEFF(J)=C(J,1) 100 CONTINUE CC-----IF DIAGNOSTICS ARE WANTED, THEN CALCULATE THE u_i IF (JDIAG.EQ.0) GOTO 200 DO 110 JNC=1,NCAS UHYB(JNC)=ABS(RESDU(JNC))/FCTW 110 CONTINUE ELSE CC-----UPDATE THE "BEST" ESTIMATE AND OBJECTIVE FUNCTION VALUE CC-----IF NECESSARY IF (FCTW.LT.COEFF(NVAD)) THEN COEFF(NVAD)=FCTW DO 140 J=1,NVAR JBEST(J)=JSUBS(J) 140 COEFF(J)=C(J,1) ENDIF IF (JDIAG.EQ.0) GOTO 200 DO 130 JNC=1,NCAS UH=ABS(RESDU(JNC))/FCTW IF(UHYB(JNC).LT.UH) UHYB(JNC)=UH 130 CONTINUE ENDIF CC-----IF FCTW IS LESS THAN PREC (A PRECISION CONSTANT) THEN ONE IS CC-----RESORTED TO THE EXACT FIT CASE 200 IF (JHYB.EQ.0) THEN COEFF(NVAD)=FCTW DO 150 J=1,NVAR JBEST(J)=JSUBS(J) 150 COEFF(J)=C(J,1) CALL RESTT(COEFF,1,0,NVAR,NCAS,NVAD,NZWE,X,Y,RESDU,WEIGHTS, 1 XMED,XMAD,MAXP1,MAXN,0,SCAL,QUAN,PREC) AVW=NZWE/AL RETURN ENDIF 10 CONTINUE 210 CONTINUE CC-----FCTW : LOWEST VALUE OF OBJECTIVE FUNCTION CC-----PRELSCAL : PRELIMINARY SCALE ESTIMATE FCKW=COEFF(NVAD) PRELSCAL=FCKW*FACTOR FINITF=-10.0/(AL-ANVAR)*(JQU/AL)+(10.0/(AL-ANVAR)+1.0) IF (MROB.EQ.1) PRELSCAL=PRELSCAL*FINITF COEFF(NVAD)=PRELSCAL CC-----COEFFICIENT OF DETERMINATION CORRESPONDING TO THE LMS/LTS RSQ=PRELSCAL DO 230 JNC=1,NCAS RESDU(JNC)=Y(JNC) 230 CONTINUE IF (JCST.NE.0) THEN CALL RANGS(RESDU,NCAS) IF (MROB.EQ.1) THEN CALL LMSLOC(RESDU,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) ELSE CALL LTSLOC(RESDU,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) ENDIF RSY=BSTD ELSE CC-----REGRESSION TROUGH ORIGIN, FIRST TRANSFORM Y DO 240 JNC=1,NCAS RESDU(JNC)=RESDU(JNC)*XMAD(NVAD)+XMED(NVAD) 240 CONTINUE RSQ=RSQ*XMAD(NVAD) IF (MROB.EQ.1) THEN DO 295 JNC=1,NCAS AW2(JNC)=ABS(RESDU(JNC)) 295 CONTINUE RSY=PULL(AW,MAXN,AW2,NCAS,JQU)*FACTOR ELSE DO 297 JNC=1,NCAS AW2(JNC)=RESDU(JNC)**2 297 CONTINUE CALL RANGS(AW2,NCAS) RSY=0 DO 298 J=1,JQU RSY=RSY+AW2(J) 298 CONTINUE RSY=SQRT(RSY/JQU)*FACTOR ENDIF ENDIF IF (MROB.EQ.1) RSY=RSY*FINITF RSQ=1.0-(RSQ/RSY)**2 IF (RSQ.LT.0.0) RSQ=0.0 IF (RSQ.GT.1.0) RSQ=1.0 CC-----COMPUTING THE FINAL SCALE ESTIMATE OF THE LMS OR THE LTS CC-----AND THE RESULTING WEIGHTS CALL RESTT(COEFF,1,0,NVAR,NCAS,NVAD,NZWE, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,MAXP1,MAXN,1,SCAL,PRELSCAL,PREC) CC-----COEFF(NVAD) : FINAL SCALE ESTIMATE COEFF(NVAD)=SCAL AVW=NZWE/AL 8010 FORMAT(I4,48h percent of the calculations have been executed.) RETURN END CC ------------------------------------------------------------------------ CC PRINTING LMS/LTS ESTIMATES AND SOME OPTIONS CC ------------------------------------------------------------------------ SUBROUTINE PRTLMS(NVAR,NCAS,JCST,NFAC,NVAD,MAXP1,MAXN,MAXP, 1 XMED,XMAD,COEFF,AW,RESDU,WEIGHTS,X,Y,JNDEX,UHYB,NZWE,RSQ,JREG, 1 JBEST,JPRT,JPLT,NREP,JERD,JLP,JQU,JHYB,JDIAG,NSTOP,JHEAD,LAB, 1 PREC,LUB,FCKW,MROB,ALGO,PRELSCAL) DIMENSION X(MAXP1,MAXN),Y(MAXN),XMED(MAXP1),XMAD(MAXP1) DIMENSION COEFF(MAXP),AW(MAXN),RESDU(MAXN) DIMENSION WEIGHTS(MAXN),UHYB(MAXN) INTEGER JNDEX(MAXN) INTEGER JBEST(MAXP) CHARACTER ALGO CHARACTER*10 LAB(MAXP1) CHARACTER*60 JHEAD CHARACTER*3 NAME NAME=' ' IF (MROB.EQ.1) THEN NAME='LMS' ELSE NAME='LTS' ENDIF WRITE(LUB,8000) IF (MROB.EQ.1) THEN WRITE(LUB,8010) WRITE(LUB,8020) JQU ELSE WRITE(LUB,8015) WRITE(LUB,8025) JQU ENDIF JBREAK=NBREAK(JQU,NCAS,NVAR) JDEFAUL=(NCAS+NVAR+1)/2 IF (JBREAK.EQ.NBREAK(JDEFAUL,NCAS,NVAR)) THEN WRITE(LUB,8028) JBREAK,NCAS,NVAR ELSE WRITE(LUB,8027) JBREAK ENDIF NSTOP=0 IF (ALGO.EQ.'R') THEN WRITE(LUB,8030) JLP,NVAR,NCAS IF (JERD.NE.0) THEN IF (JERD.EQ.JLP) THEN WRITE(LUB,8035) NSTOP=1 RETURN ELSE WRITE(LUB,8040) JLP,JERD ENDIF ELSE WRITE(LUB,8045) JLP ENDIF ELSE IF (JHYB.NE.0) THEN WRITE(LUB,8050) NREP,NVAR,NCAS ELSE WRITE(LUB,8055) JLP,NVAR,NCAS ENDIF IF (JERD.NE.0) THEN IF (JERD.EQ.NREP) THEN WRITE(LUB,8035) NSTOP=1 RETURN ELSE IF (JHYB.NE.0) THEN WRITE(LUB,8040) NREP,JERD ELSE WRITE(LUB,8040) JLP,JERD ENDIF ENDIF ELSE IF (JHYB.NE.0) WRITE(LUB,8045) NREP ENDIF ENDIF WRITE(LUB,8065) (JBEST(J),J=1,NVAR) CC-----LMS/LTS IN THE CASE OF EXACT FIT IF (JHYB.EQ.0) THEN WRITE(LUB,8070) DO 10 J=1,NFAC 10 WRITE(LUB,8080) LAB(J),COEFF(J) IF (JCST.EQ.0) WRITE(LUB,8080) LAB(NVAR),COEFF(NVAR) IF (JCST.NE.0) WRITE(LUB,8090) COEFF(NVAR) WRITE(LUB,8100) COEFF(NVAD) WRITE(LUB,8110) NZWE JREG=2 CC-----RESIDUAL PLOT FOR THE LMS/LTS IN THE CASE OF EXACT FIT CALL RDUAL(COEFF,0,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC, 1 JREG, X,Y,RESDU,WEIGHTS,XMED,XMAD,JPLT,AW,JNDEX, 1 MAXP1,MAXN,JHEAD,LAB,MROB) NSTOP=1 RETURN ENDIF CC-----LMS/LTS IN THE GENERAL CASE WRITE(LUB,8130) DO 20 J=1,NFAC 20 WRITE(LUB,8080) LAB(J),COEFF(J) IF (JCST.EQ.0) WRITE(LUB,8080) LAB(NVAR),COEFF(NVAR) IF (JCST.NE.0) WRITE(LUB,8090) COEFF(NVAR) WRITE(LUB,8190) NAME,FCKW WRITE(LUB,8200) NAME,PRELSCAL WRITE(LUB,8210) RSQ WRITE(LUB,8215) NAME,COEFF(NVAD) WRITE(LUB,8235) WRITE(LUB,8240) NCAS,JQU,NVAR IF (MROB.EQ.1) THEN WRITE(LUB,8250) WRITE(LUB,8300) ELSE WRITE(LUB,8260) WRITE(LUB,8310) ENDIF IF (JCST.NE.0) THEN WRITE(LUB,8320) ELSE WRITE(LUB,8325) ENDIF WRITE(LUB,8335) NAME WRITE(LUB,8235) JREG=2 CC-----RESIDUAL PLOTS IF REQUESTED CALL RDUAL(COEFF,1,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,JPLT,AW,JNDEX,MAXP1,MAXN, 1 JHEAD,LAB,MROB) CC-----RESISTANT DIAGNOSTICS IF (JHYB.NE.0.AND.JDIAG.EQ.1) THEN WRITE(LUB,8350) WRITE(LUB,8360) NAME UH=AMDAN(AW,MAXN,UHYB,NCAS) DO 40 JNC=1,NCAS RDI=UHYB(JNC)/UH WRITE(LUB,8390) JNDEX(JNC),UHYB(JNC),RDI 40 CONTINUE IF (MROB.EQ.1) THEN WRITE(LUB,8400) ELSE WRITE(LUB,8410) ENDIF WRITE(LUB,8420) ENDIF IF(JCST.EQ.1) THEN WRITE(LUB,8430) NAME WRITE(LUB,8440) ENDIF 8000 FORMAT(/1X,78(1H*)/) 8010 FORMAT(' LEAST MEDIAN OF SQUARES REGRESSION'/, 1 1X,34(1H*)/) 8015 FORMAT(' LEAST TRIMMED SQUARES REGRESSION'/, 1 1X,32(1H*)/) 8020 FORMAT(' The estimator minimizes the ',I4, 1 'th ordered squared residual,') 8027 FORMAT(' hence its breakdown value is ',I2,'%.') 8028 FORMAT(' hence its breakdown value is ',I2,'% (= highest', 1 ' possible value '/' for n = ',I4,' cases and p = ',I3, 1 ' coefficients).') 8025 FORMAT(' The estimator minimizes the sum of the ',I4, 1 ' smallest squared residuals,') 8030 FORMAT(/' The program has drawn ',I8,' random subsets of ',I3, 1 ' cases out of ',I5,'.') 8035 FORMAT(//' Each of these subsets', 1 ' gave a singular system of equations.'/ 1 48h This means that there is an enormous problem of, 1 30h collinearity in the data set./ 1 45h Try the analysis again in a lower dimension./) 8040 FORMAT(' Among these ',I8,' subsets,',I8,' were discarded'/ 1 ' because they gave a singular system of equations.') 8045 FORMAT(' Each of these ',I8,' subsets gave a nonsingular system', 1 ' of equations.') 8050 FORMAT(/' The program has considered all ',I7,' subsets of ',I3, 1 ' cases out of ',I5,'.') 8055 FORMAT(/' The program has drawn ',I8,' subsets of ',I3, 1 ' cases out of ',I5,'.') 8065 FORMAT(/' The best subset consisted of the observations'/14I5//) 8070 FORMAT(/' More than the half of the data lie on the same ', 1 'hyperplane,'/' with coefficients = '//) 8080 FORMAT(3X,A10,F20.5) 8090 FORMAT(3X,10h constant,F20.5) 8100 FORMAT(/15h Scale estimate,10X,2H= ,F15.5/) 8110 FORMAT(10h There are,2X,I4,' points on the hyperplane.') 8130 FORMAT(//5X,8hvariable,9X,11hcoefficient/3X,30(1H-)) 8190 FORMAT(//1X,A3,' objective function = ',F15.5) 8200 FORMAT(/' Preliminary ',A3,' scale = ',F15.5) 8210 FORMAT(/' Robust R-squared (coefficient of determination) = ', 1 F7.5) 8215 FORMAT(/' Final ',A3,' scale estimate = ',F15.5) 8235 FORMAT(/' ---------------------------------------------------', 1 '-------------------------') 8240 FORMAT(/' Corresponding formulas: (for this data set ', 1 'n = ',I5,', h = ',I5,', p = ',I3,')') 8250 FORMAT(// 1 ' LMS objective = obj(x ,y ) = minimum |y - x *theta|'/ 1 ' i i theta i i h:n') 8260 FORMAT(// 1 ' | 1 h ', 1 ' 2 |'/ 1 ' LTS objective = obj(x ,y ) = minimum squareroot | - sum', 1 ' |y - x *theta| |'/ 1 ' i i theta | h j=1', 1 ' i i j:n |') 8300 FORMAT(// 1 ' Preliminary LMS scale = s(x ,y ) = constant*obj(x ,y )'/ 1 ' i i i i ') 8310 FORMAT(// 1 ' Preliminary LTS scale = s(x ,y ) = constant*obj(x ,y )'/ 1 ' i i i i ') 8320 FORMAT(// 1 ' | s(x ,y ) | 2'/ 1 ' | i i | '/ 1 ' Robust R-squared = 1 - | ---------- | '/ 1 ' | s(1,y ) | '/ 1 ' | i | ') 8325 FORMAT(// 1 ' | s(x ,y ) | 2'/ 1 ' | i i | '/ 1 ' Robust R-squared = 1 - | ---------- | '/ 1 ' | s(0,y ) | '/ 1 ' | i | ') 8335 FORMAT(// 1 ' | n 2 |'/ 1 ' | sum w *r |'/ 1 ' | i=1 i i |'/ 1 ' Final ',A3,' scale estimate = squareroot | -------------- |'/ 1 ' | n |'/ 1 ' | (sum w ) - p |'/ 1 ' | i=1 i |') 8350 FORMAT(/' ---------------------------------------------------', 1 '-------------------------'//' Resistant diagnostic: '/ 1 ' *********************'/) 8360 FORMAT(' This diagnostic is based on the ',A3,' and becomes', 1 ' larger than 2.5 (say) '/' in regression outliers as well', 1 ' as in good leverage points'/' and in bad leverage points.'// 1 1X,' ',3X,13h u ,3X,10hresistant / 1 1X,' i ',12X,'i',6X,10hdiagnostic/1X,34(1H-)) 8390 FORMAT(1X,I4,3X,F13.4,3X,F10.4) 8400 FORMAT(// 1 ' | y - x *theta | '/ 1 ' i i '/ 1 ' where u = maximum ------------------- '/ 1 ' i theta | y - x *theta | '/ 1 ' j j h:n') 8410 FORMAT(// 1 ' | y - x *theta | '/ 1 ' i i '/ 1 ' where u = maximum -------------------------------------- '/ 1 ' i theta | 1 h 2 |'/ 1 ' squareroot | - sum |y - x *theta| |'/ 1 ' | h j=1 k k j:n |') 8420 FORMAT(// 1 ' n '/ 1 ' and resistant diagnostic = u / median (u ) '/ 1 ' i i j=1 j ') 8430 FORMAT(//' ---------------------------------------------------', 1 '-------------------------'// 1 ' A complete diagnosis of the cases can be obtained as', 1 ' follows:'// 1 ' (1) The standardized ',A3,' residuals provided by this', 1 ' program'/ 1 ' indicate regression outliers, i.e. cases for which'// 1 ' |res /s| > 2.5 (roughly).'/ 1 ' i '// 1 ' (2) Run the program MINVOL on the data set formed by the', 1 ' explanatory '/ 1 ' variables only. This yields the minimum', 1 ' volume ellipsoid (MVE)'/ 1 ' scatter matrix and the corresponding robust distances RD '/ 1 ' with their cutoff value c. i'/ 1 ' Leverage points (= x-outliers) are then cases with'// 1 ' RD > c (roughly).'/ 1 ' i ') 8440 FORMAT(/' (3) Combining this information classifies the cases', 1 ' into four types:'// 1 ' (a) regular observation: RD <= c and |res /s| <= 2.5'/ 1 ' i i '// 1 ' (b) vertical outlier : RD <= c and |res /s| > 2.5'/ 1 ' i i '// 1 ' (c) good leverage point: RD > c and |res /s| <= 2.5'/ 1 ' i i '// 1 ' (d) bad leverage point : RD > c and |res /s| > 2.5'/ 1 ' i i '// 1 ' For more background and some examples see'// 1 ' Rousseeuw, P.J., and van Zomeren, B.C. (1990),'/ 1 ' Unmasking Multivariate Outliers and Leverage Points,'/ 1 ' Journal of the American Statistical Association, 85,', 1 ' 633-639'// 1 ' where a plot of res /s versus RD is proposed.'/ 1 ' i i ') RETURN END CC ------------------------------------------------------------------------ CC *RDAT*: READ DATA AND OPTIONS CC ------------------------------------------------------------------------ SUBROUTINE RDAT(NCAS,NVAR,JCST,JPRT,NFAC,NVAD,VALMS, 1 X,Y,JPLT,JMISS,PREC,JPLACE,LAB,XMED,XMAD,AW,AW2,JNDEX, 1 RESDU,WEIGHTS,MAXP1,MAXN,MVAL,MAXP,LUA,LUB,LUC,JHEAD, 1 YNSAVE,JDIAG,NSTOP,FNAMEA,FNAMEB,FNAMEC,MROB,FACLMS,FACLTS) DIMENSION X(MAXP1,MAXN),Y(MAXN) DIMENSION AW(MAXN),AW2(MAXN),RESDU(MAXN) DIMENSION WEIGHTS(MAXN),XMED(MAXP1) DIMENSION XMAD(MAXP1),VALMS(MAXP1) DIMENSION FACLMS(11),FACLTS(11) INTEGER JMISS(MAXP1),JNDEX(MAXN),JPLACE(MAXP1) CHARACTER YN,YNOK,YNSAVE,YNM CHARACTER*30 FNAMEA,FNAMEB,FNAMEC CHARACTER*10 LAB(MAXP1) CHARACTER*60 JFMT,JHEAD CHARACTER*10 LABJ CHARACTER HELP(10) LOGICAL EXISTA EQUIVALENCE (LABJ,HELP) JDIAG=0 LOCA=0 WRITE (*,8000) CC-----GIVE THE NUMBER OF CASES 10 WRITE (*,8010) 11 READ(*,*,ERR=11)NCAS IF (NCAS.GT.MAXN.OR.NCAS.LE.1) THEN IF (NCAS.GT.MAXN) WRITE (*,8020) MAXN IF (NCAS.LE.1) WRITE(*,8030) GOTO 10 ENDIF CC-----CONSTANT TERM OR NOT WRITE (*,8040) 20 READ (*,9000) YN IF (YN.EQ.'y') YN='Y' IF (YN.EQ.'n') YN='N' IF (YN.NE.'Y'.AND.YN.NE.'N') THEN WRITE (*,8050) GOTO 20 ELSE IF (YN.EQ.'Y') JCST=1 IF (YN.EQ.'N') JCST=0 ENDIF CC-----TOTAL NUMBER OF VARIABLES IN THE DATA SET WRITE(*,8060) 30 READ(*,*,ERR=30)JVARS IF (JVARS.LT.1.OR.JVARS.GT.50) THEN WRITE(*,8050) GOTO 30 ENDIF IF (JVARS.EQ.1) THEN CC-----THE PROBLEM IS REDUCED TO LOCATION LOCA=1 NVAR=1 JCST=1 WRITE(*,8070) JPLACE(MAXP1)=1 ELSE CC-----WHICH VARIABLE WILL BE TAKEN AS THE RESPONSE VARIABLE WRITE(*,8080) JVARS 40 READ(*,*,ERR=40)JPLACE(MAXP1) IF (JPLACE(MAXP1).LT.1.OR.JPLACE(MAXP1).GT.JVARS) THEN WRITE(*,8050) GOTO 40 ENDIF ENDIF CC-----GIVE A LABEL FOR THIS VARIABLE WRITE(*,8090) READ(*,9010) LABJ CALL MOVE(HELP) LAB(MAXP1)=LABJ IF (LOCA.EQ.1) GOTO 100 IF (JCST.EQ.0) THEN CC-----JWR = MAXIMAL NUMBER OF EXPLANATORY VARIABLES POSSIBLE JWR=MAXP IF((JVARS-1).LT.MAXP) JWR=JVARS-1 ELSE JWR=MAXP-1 IF ((JVARS-1).LT.MAXP) JWR=JVARS-1 ENDIF CC-----HOW MANY EXPLANATORY VARIABLES DO YOU WANT TO USE WRITE(*,8100) JWR 60 READ(*,*,ERR=60)NVAR IF (NVAR.EQ.0) THEN WRITE(*,8070) LOCA=1 NVAR=1 JCST=1 GOTO 100 ENDIF IF (NVAR.LT.1.OR.NVAR.GT.JWR) THEN WRITE(*,8050) GOTO 60 ENDIF IF (NVAR.NE.(JVARS-1)) GOTO 80 IF (JCST.EQ.1) NVAR=JVARS IF (JCST.EQ.0) NVAR=JVARS-1 IF (NVAR.GT.MAXP) WRITE (*,8115) MAXP IF (NCAS.GT.MAXN.OR.NVAR.GT.MAXP) GOTO 10 IF (NCAS.LE.(NVAR*2)) THEN CC----- TOO MANY CASES WRITE (*,8110) NCAS,NVAR IF (JCST.EQ.1) WRITE(*,8120) GOTO 10 ENDIF NVAD=NVAR+1 JPLACE(NVAD)=JPLACE(MAXP1) LAB(NVAD)=LAB(MAXP1) JP=0 WRITE(*,8130) DO 70 J=1,NVAR JP=JP+1 IF (.NOT.(JCST.EQ.1.AND.J.EQ.NVAR)) THEN IF (J.EQ.JPLACE(MAXP1)) JP=JP+1 JPLACE(J)=JP WRITE(*,8140) J,JP READ(*,9010) LABJ CALL MOVE(HELP) LAB(J)=LABJ ENDIF 70 CONTINUE GOTO 135 80 IF (JCST.EQ.1) NVAR=NVAR+1 IF (NCAS.LE.(NVAR*2)) THEN WRITE (*,8110) NCAS,NVAR IF (JCST.EQ.1) WRITE(*,8120) GOTO 10 ENDIF 100 NVAD=NVAR+1 JPLACE(NVAD)=JPLACE(MAXP1) LAB(NVAD)=LAB(MAXP1) IF (LOCA.EQ.1) GOTO 135 CC-----EXPLANATORY VARIABLES WRITE(*,8130) DO 110 J=1,NVAR JTT=J-1 IF (JCST.EQ.1.AND.J.EQ.NVAR) GOTO 110 120 WRITE(*,8150) J READ(*,9020) JPLACE(J),LABJ CALL MOVE(HELP) LAB(J)=LABJ IF (JPLACE(J).LT.1.OR.JPLACE(J).GT.JVARS) GOTO 120 IF (JPLACE(J).EQ.JPLACE(NVAD)) THEN WRITE (*,8160) GOTO 120 ENDIF IF (J.EQ.1) GOTO 110 DO 130 JK=1,JTT IF (JPLACE(JK).NE.JPLACE(J))GOTO 130 WRITE(*,8160) GOTO 120 130 CONTINUE 110 CONTINUE CC-----LMS OR LTS? 135 WRITE (*,8165) 136 READ (*,9000) YN IF(.NOT.(YN.EQ.'1'.OR.YN.EQ.'2')) THEN WRITE (*,8050) GOTO 136 ENDIF MROB=JNTGR(YN) CC-----HOW MUCH OUTPUT? IF (LOCA.NE.1) THEN WRITE (*,8170) 150 READ (*,9000) YN IF(.NOT.(YN.EQ.'0'.OR.YN.EQ.'1'.OR.YN.EQ.'2')) THEN WRITE (*,8050) GOTO 150 ENDIF ELSE WRITE (*,8175) 155 READ (*,9000) YN IF(.NOT.(YN.EQ.'0'.OR.YN.EQ.'1')) THEN WRITE (*,8050) GOTO 155 ENDIF ENDIF JPRT=JNTGR(YN) CC-----DO YOU WANT TO LOOK AT THE RESIDUALS? IF (LOCA.NE.1) THEN WRITE(*,8180) ELSE WRITE (*,8190) ENDIF 160 READ (*,9000) YN IF(.NOT.(YN.EQ.'0'.OR.YN.EQ.'1'.OR.(YN.EQ.'2'.AND.LOCA.NE.1) 1 .OR.(YN.EQ.'3'.AND.LOCA.NE.1) ) ) THEN WRITE (*,8050) GOTO 160 ENDIF JPLT=JNTGR(YN) IF (JCST.EQ.1.AND.NVAR.EQ.1.AND.JPLT.NE.0) JPLT=2 CC-----DO YOU WANT TO COMPUTE OUTLIER DIAGNOSTICS IF (LOCA.NE.1) THEN WRITE(*,8200) 162 READ (*,9000) YN IF (YN.EQ.'y') YN='Y' IF (YN.EQ.'n') YN='N' IF (YN.NE.'Y'.AND.YN.NE.'N') THEN WRITE (*,8050) GOTO 162 ENDIF IF (YN.EQ.'Y') THEN JDIAG=1 ELSE JDIAG=0 ENDIF ENDIF CC-----NAME OF THE FILE CONTAINING THE DATA 165 WRITE(*,8210) READ(*,9030) FNAMEA IF (FNAMEA.EQ.'Key'.OR.FNAMEA.EQ.'key'.OR.FNAMEA.EQ.'KEY') THEN FNAMEA='CON' ELSE INQUIRE(FILE=FNAMEA,EXIST=EXISTA,IOSTAT=NER) IF (NER.NE.0.OR.(.NOT.EXISTA)) THEN IF (NER.EQ.7.OR.NER.EQ.10.OR.NER.EQ.1027.OR.NER.EQ.1032.OR. 1 NER.EQ.1033.OR.NER.EQ.1002.OR.NER.EQ.1003.OR. 1 NER.EQ.1004.OR.NER.EQ.1005.OR.NER.EQ.1007.OR. 1 (.NOT.EXISTA)) THEN WRITE(*,8220) GOTO 165 ELSE WRITE(*,8230) NER NSTOP=1 RETURN ENDIF ENDIF ENDIF IF (FNAMEA.EQ.'CON') THEN WRITE(*,8240) CC-----DO YOU WANT TO SAVE THE DATA 170 READ(*,9000) YNSAVE IF (YNSAVE.EQ.'y') YNSAVE='Y' IF (YNSAVE.EQ.'n') YNSAVE='N' IF (YNSAVE.NE.'Y'.AND.YNSAVE.NE.'N') THEN WRITE (*,8050) GOTO 170 ENDIF IF (YNSAVE.EQ.'Y') THEN 175 WRITE(*,8250) READ(*,9030) FNAMEC INQUIRE(FILE=FNAMEC,EXIST=EXISTA) IF (EXISTA) THEN WRITE(*,8225) GOTO 175 ENDIF OPEN(LUC,FILE=FNAMEC,STATUS='NEW') ENDIF ENDIF OPEN(LUA,FILE=FNAMEA) CC-----WHERE DO YOU WANT THE OUTPUT 185 WRITE(*,8260) READ(*,9030) FNAMEB IF (FNAMEB.EQ.'con'.OR.FNAMEB.EQ.'Con') FNAMEB='CON' IF (FNAMEB.EQ.'prn'.OR.FNAMEB.EQ.'Prn') FNAMEB='PRN' IF (FNAMEB.EQ.'CON'.OR.FNAMEB.EQ.'PRN') THEN OPEN(LUB,FILE=FNAMEB) ELSE INQUIRE(FILE=FNAMEB,EXIST=EXISTA) IF (EXISTA) THEN WRITE(*,8225) GOTO 185 ENDIF OPEN(LUB,FILE=FNAMEB,STATUS='NEW') ENDIF CC-----PLEASE ENTER A TITLE FOR THE OUTPUT WRITE (*,8270) READ (*,9040) JHEAD NFAC=NVAR-1 CC-----DO YOU WANT TO READ THE DATA IN FREE FORMAT WRITE (*,8280) 180 READ (*,9000) YN IF (YN.EQ.'y') YN='Y' IF (YN.EQ.'n') YN='N' IF (YN.NE.'Y'.AND.YN.NE.'N') THEN WRITE (*,8050) GOTO 180 ENDIF IF (YN.NE.'Y') THEN CC-----READ FORMAT FOR THE INPUT WRITE(*,8290) READ (*,9040) JFMT ENDIF CC-----OPTION FOR THE TREATMENT OF MISSING VALUE IF (LOCA.NE.1) THEN WRITE (*,8310) ELSE WRITE (*,8320) ENDIF 200 READ (*,9000) YNOK IF (.NOT.(YNOK.EQ.'0'.OR.YNOK.EQ.'1'.OR. 1 (YNOK.EQ.'2'.AND.LOCA.NE.1) ) ) THEN WRITE (*,8050) GOTO 200 ENDIF MVAL=JNTGR(YNOK) CC-----ABSTRACT OF ALL THE OPTIONS WRITE(*,8325) WRITE(*,8326) JHEAD WRITE(*,8330) NCAS IF (LOCA.EQ.1) THEN WRITE(*,8340) LAB(NVAD) WRITE(*,8345) ELSE IF (JCST.EQ.0) WRITE(*,8335) NVAR IF (JCST.EQ.1) WRITE(*,8335) NFAC WRITE(*,8340) LAB(NVAD) IF (JCST.EQ.0) WRITE(*,8355) IF (JCST.EQ.1) WRITE(*,8350) ENDIF IF (MROB.EQ.1) WRITE(*,8360) IF (MROB.EQ.2) WRITE(*,8365) IF (FNAMEA.NE.'CON') WRITE(*,8370) FNAMEA IF (FNAMEA.EQ.'CON') WRITE(*,8380) IF (YNSAVE.EQ.'Y') WRITE(*,8390) FNAMEC IF (YN.EQ.'Y') WRITE(*,8410) IF (YN.EQ.'N') WRITE(*,8420) JFMT IF (JPRT.EQ.0) WRITE (*,8430) IF (JPRT.EQ.1.AND.LOCA.NE.1) WRITE (*,8440) IF (JPRT.EQ.2.OR.(JPRT.EQ.1.AND.LOCA.EQ.1)) WRITE (*,8450) IF (JPLT.EQ.0) WRITE(*,8460) IF (JPLT.EQ.1) WRITE(*,8470) IF (JPLT.EQ.2) WRITE(*,8480) IF (JPLT.EQ.3) WRITE(*,8490) IF (MVAL.EQ.0) WRITE(*,8520) IF (MVAL.EQ.1.AND.LOCA.NE.1) WRITE(*,8530) IF (MVAL.EQ.1.AND.LOCA.EQ.1) WRITE(*,8535) IF (MVAL.EQ.2) WRITE(*,8540) WRITE(*,8550) FNAMEB WRITE (*,8560) 210 READ (*,9000) YNOK IF (YNOK.EQ.'y') YNOK='Y' IF (YNOK.EQ.'n') YNOK='N' IF (YNOK.NE.'Y'.AND.YNOK.NE.'N') THEN WRITE (*,8050) GOTO 210 ENDIF IF (YNOK.NE.'Y') GOTO 10 IF (LOCA.EQ.1) 1 CALL LCAT(NCAS,NVAR,JCST,JPRT,NVAD,X,Y,RESDU,WEIGHTS,PREC,MROB, 1 XMED,XMAD,NZWE,AVW,JPLT,AW,JNDEX,MAXP1,MAXN,MVAL,LUA,LUB,LUC, 1 JREG,JHEAD,FNAMEA,FNAMEB,FNAMEC,YNSAVE,LAB,JFMT,JVARS,YN,JPLACE, 1 AW2,FACLMS,FACLTS) IF (MVAL.NE.0) THEN DO 220 J=1,NVAD IF (J.EQ.1) THEN WRITE (*,8570) 225 READ (*,9000) YNM IF (YNM.EQ.'y') YNM='Y' IF (YNM.EQ.'n') YNM='N' IF (YNM.NE.'Y'.AND.YNM.NE.'N') THEN WRITE(*,8050) GOTO 225 ENDIF IF (YNM.EQ.'Y'.OR.YNM.EQ.'y') THEN YNM='Y' WRITE (*,8580) 226 READ(*,*,ERR=226)CODE ENDIF ENDIF IF (YNM.NE.'Y') THEN IF (J.EQ.NVAR.AND.JCST.EQ.1) GOTO 220 IF (J.NE.NVAD) WRITE (*,8590) LAB(J) IF (J.EQ.NVAD) WRITE (*,8600) 230 READ (*,9000) YNOK IF (YNOK.EQ.'y') YNOK='Y' IF (YNOK.EQ.'n') YNOK='N' IF (YNOK.NE.'Y'.AND.YNOK.NE.'N') THEN WRITE (*,8050) GOTO 230 ENDIF ELSE YNOK='Y' ENDIF IF (YNOK.NE.'Y') THEN JMISS(J)=0 ELSE JMISS(J)=1 CC-----READ THE VALUE WHICH HAS TO BE INTERPRETED AS THE MISSING VALUE CODE IF (YNM.NE.'Y') THEN WRITE(*,8610) 231 READ(*,*,ERR=231)VALMS(J) ELSE VALMS(J)=CODE ENDIF ENDIF 220 CONTINUE ENDIF CC-----ENTER THE DATA FOR EACH CASE IF (FNAMEA.EQ.'CON') WRITE(*,8615) DO 250 JNC=1,NCAS IF (FNAMEA.EQ.'CON') THEN IF (JNC.EQ.1) WRITE(*,8620) JNC IF (JNC.NE.1) WRITE(*,8630) JNC ENDIF JNDEX(JNC)=JNC IF (JCST.EQ.0) THEN IF (YN.EQ.'Y') THEN READ(LUA,*)(AW(J),J=1,JVARS) DO 260 J=1,NVAD JH=JPLACE(J) X(J,JNC)=AW(JH) 260 CONTINUE IF (YNSAVE.EQ.'Y') WRITE(LUC,*)(AW(J),J=1,JVARS) ELSE READ(LUA,JFMT)(AW(J),J=1,JVARS) DO 270 J=1,NVAD JH=JPLACE(J) X(J,JNC)=AW(JH) 270 CONTINUE IF (YNSAVE.EQ.'Y') WRITE(LUC,*)(AW(J),J=1,JVARS) ENDIF ELSE X(NVAR,JNC)=1.0 IF (YN.EQ.'Y') THEN READ(LUA,*)(AW(J),J=1,JVARS) DO 280 J=1,NVAD IF (.NOT.(J.EQ.NVAR.AND.JCST.EQ.1)) THEN JH=JPLACE(J) X(J,JNC)=AW(JH) ENDIF 280 CONTINUE IF (YNSAVE.EQ.'Y') WRITE(LUC,*)(AW(J),J=1,JVARS) ELSE READ(LUA,JFMT)(AW(J),J=1,JVARS) DO 290 J=1,NVAD IF (.NOT.(J.EQ.NVAR.AND.JCST.EQ.1)) THEN JH=JPLACE(J) X(J,JNC)=AW(JH) ENDIF 290 CONTINUE IF (YNSAVE.EQ.'Y') WRITE(LUC,*) (AW(J),J=1,JVARS) ENDIF ENDIF 250 CONTINUE IF (YNSAVE.EQ.'Y') CLOSE(LUC,STATUS='KEEP') IF (YNSAVE.EQ.'Y') WRITE(*,8395) FNAMEC CC-----FORMATS 8000 FORMAT(/////////////////30X,19(1H*)/ 1 30X,19H* P R O G R E S S */30X,19(1H*)///////) 8010 FORMAT(/36h Enter the number of cases please : ,$) 8020 FORMAT(/34h There are too many cases (at most,I5,2H ), 1 24h according to the limits 1 /58h of the arrays in the program. The user has to adapt these 1 /59h limits such that the concerning limits are greater than or 1 /48h equal to the number of cases in his data set or 1 /37h another regression can be performed.//) 8030 FORMAT(/27h There are not enough cases, 1 27h (at least 2 are required).//) 8040 FORMAT(/47h Do you want a constant term in the regression?/ 1 20h Answer yes or no : ,$) 8050 FORMAT(42h Not allowed ! Enter your choice again : ,$) 8060 FORMAT(/46h What is the total number of variables in your, 1 10h data set?/1X,55(1H-)/ 1 41h Please give a number between 1 and 50 : ,$) 8070 FORMAT(//' The problem is reduced to estimating univariate', 1 ' location and scale.') 8080 FORMAT(/51h Which variable do you choose as response variable?/ 1 1X,50(1H-)/14h Out of these ,I4,21h give its position : ,$) 8090 FORMAT(/31h Give a label for this variable, 1 27h (at most 10 characters) : ,$) 8100 FORMAT(/50h How many explanatory variables do you want to use, 1 17h in the analysis?/1X,66(1H-)/ 1 10h (at most ,I4,5H ) : ,$) 8110 FORMAT(//46h Too many coefficients according to the number, 1 10h of cases.//16h Number of cases,8X,2H= ,I5/ 1 26h Number of coefficients = ,I5/25h The number of cases must, 1 37h be twice the number of coefficients.) 8115 FORMAT(//41h There are too many coefficients (at most,I3,2H ), 1 24h according to the limits/30h of the arrays in the program., 1 28h The user has to adapt these/21h limits such that the, 1 38h concerning limits are greater than or/ 1 50h equal to the number of coefficients in the model.//) 8120 FORMAT(31h (Including the constant term!)) 8130 FORMAT(//40h Explanatory variables : position, 1 32h label (at most 10 characters)/ 1 1X,32(1h-),4(1h+),6(1h-),10(1h+),19(1h-)) 8140 FORMAT(9h Number ,I4,15X,1H:,4X,I4,5X,' ',$) 8150 FORMAT(9h Number ,I4,15X,': ',$) 8160 FORMAT(/51h This position has already been choosen for another, 1 10h variable./35h Enter the right position please : ) 8165 FORMAT(/' Which robust method do you want to use? '/ 1 ' --------------------------------------- '/ 1 ' 1 = Least median of squares (LMS)'/ 1 ' 2 = Least trimmed squares (LTS)'/ 1 ' Enter your choice : ',$) 8170 FORMAT(/29h How much output do you want?/1X,28(1H-)/ 1 17h 0 = small output,7X,26h: limited to basic results/ 1 25h 1 = medium-sized output:, 1 53h also includes a table with the observed values of y,/ 1 25X,50h the estimates of y, the residuals and the weights/ 1 17h 2 = large output,7X, 1 45h: also includes the raw and standardized data/ 1 21h Enter your choice : ,$) 8175 FORMAT(/29h How much output do you want?/1X,28(1H-)/ 1 ' 0 = small output: limited to basic results'/ 1 ' 1 = large output:', 1 ' also includes a table with the observed values of y,'/ 1 18X,' the residuals and the weights'/ 1 21h Enter your choice : ,$) 8180 FORMAT(/38h Do you want to look at the residuals?/1X,37(1H-)/ 1 22h 0 = no residual plots/ 1 56h 1 = plot of standardized residuals versus the estimated, 1 11h value of y/39h 2 = plot of the standardized residuals, 1 36h versus the index of the observation/13h 3 = performs, 1 29h both types of residual plots/21h Enter your choice : ,$) 8190 FORMAT(/38h Do you want to look at the residuals?/1X,37(1H-)/ 1 22h 0 = no residual plots/ 1 39h 1 = plot of the standardized residuals, 1 36h versus the index of the observation/ 1 21h Enter your choice : ,$) 8200 FORMAT(/' Do you want the corresponding resistant diagnostic?'/ 1 20h Answer yes or no : ,$) 8210 FORMAT(/46h Give the name of the file containing the data, 1 29h (e.g. type a:example.dat ),/8h or type, 1 51h key if you prefer to enter the data by keyboard./ 1 22h What do you choose ? ,$) 8220 FORMAT(//52h This file does not exist, please enter another one.) 8225 FORMAT(//52h This file already exists, please enter another one.) 8230 FORMAT(/22h Fortran error code : ,I8/) 8240 FORMAT(/42h Do you want to save your data in a file ?/ 1 20h Answer yes or no : ,$) 8250 FORMAT(/45h In which file do you want to save your data?/ 1 26h type e.g. b:save.dat : ,$) 8260 FORMAT(/31h Where do you want your output?/1X,30(1H-)/ 1 43h type con if you want it on the screen/ 1 44h or type prn if you want it on the printer/ 1 48h or type the name of a file (e.g. b:example.out)/ 1 23h What do you choose ? ,$) 8270 FORMAT(/36h Please enter a title for the output, 1 25h (at most 60 characters)./1X,60(1H-)/1X,$) 8280 FORMAT(/45h Do you want to read the data in free format?/ 1 1X,44(1H-)/49h This means that you only have to insert blank(s), 1 17h between numbers./25h (We advise users without, 1 45h knowledge of fortran formats to answer yes.)/ 1 28h Make your choice (yes/no): ,$) 8290 FORMAT(/33h Your desired fortran format is :/ 1 22h (between brackets and, 1 43h at most 60 characters, e.g. (2F3.0,F1.0) )) 8310 FORMAT(/53h Choose an option for the treatment of missing values/ 1 1X,52(1H-)/44h 0 = there are no missing values in the data/ 1 52h 1 = elimination of the cases for which at least one, 1 20h variable is missing/32h 2 = estimates are filled in for, 1 18h unobserved values/21h Enter your choice : ,$) 8320 FORMAT(/53h Choose an option for the treatment of missing values/ 1 1X,52(1H-)/44h 0 = there are no missing values in the data/ 1 43h 1 = elimination of the cases for which the, 1 20h variable is missing/21h Enter your choice : ,$) 8325 FORMAT(//////////35X,11('*')/35X,'* Summary *'/35X,11('*')//) 8326 FORMAT(' Title: ',A60) 8330 FORMAT(' Number of cases in the data set: ',18X,I5) 8335 FORMAT(' Number of explanatory variables in the regression: ',I5) 8340 FORMAT(' The response variable is: ',A10) 8345 FORMAT(' You have chosen univariate location and scale ', 1 'estimation.') 8350 FORMAT(' You have chosen a model with a constant term ', 1 '(regression with intercept).') 8355 FORMAT(' You have chosen a model without constant term ', 1 '(no intercept).') 8360 FORMAT(' You have chosen the least median of squares (LMS)', 1 ' method.') 8365 FORMAT(' You have chosen the least trimmed squares (LTS)', 1 ' method.') 8370 FORMAT(' Your data reside in the file: ',A30) 8380 FORMAT(41h The data will be read from the keyboard.) 8390 FORMAT(' The data will be saved in the file: ',A30) 8395 FORMAT(//' The data will be saved in the file: ',A30//) 8410 FORMAT(38h The data will be read in free format.) 8420 FORMAT(' Data input format:'/5X,A60) 8430 FORMAT(' You have chosen small output.') 8440 FORMAT(' You have chosen medium-sized output.') 8450 FORMAT(' You have chosen large output.') 8460 FORMAT(' You have chosen not to plot residuals.') 8470 FORMAT(' You have chosen a plot of residuals versus ', 1 'estimated y.') 8480 FORMAT(' You have chosen a plot of residuals versus ', 1 'their index.') 8490 FORMAT(' You have chosen a plot of residuals versus their index', 1 ' and '/ 1 ' versus estimated y.') 8520 FORMAT(' The data is assumed not to have missing values.') 8530 FORMAT(' Treatment of missing values in option 1:'/ 1 ' cases with at least one missing value will be deleted.') 8535 FORMAT(' Any case with a missing value will be deleted.') 8540 FORMAT(41h Treatment of missing values in option 2:/ 1 ' First, any case with a missing value for the response', 1 ' variable'/ 1 ' or for all explanatory variables will be deleted.'/ 1 ' For the remaining cases, each missing value is replaced by'/ 1 ' the median of the non-missing values.') 8550 FORMAT(' The output will be written on the file: ',A30) 8560 FORMAT(/27h Are all these options ok ?,13H yes or no : ,$) 8570 FORMAT(/39h Is there a unique value which is to be, 1 12h interpreted/43h as a missing measurement for any variable?/ 1 20h Answer yes or no : ,$) 8580 FORMAT(/27h Please enter this value : ,$) 8590 FORMAT(15h Does variable ,A10, 1 26h contain missing value(s)?/20h Answer yes or no : ,$) 8600 FORMAT(53h Does the response variable contain missing value(s)?/ 1 20h Answer yes or no : ,$) 8610 FORMAT(33h Enter the value of this variable, 1 31h which has to be interpreted as/ 1 26h the missing value code : ,$) 8615 FORMAT(//31h Enter your data for each case.//) 8620 FORMAT(1X,26h The data for case number ,I4,3H : ,$) 8630 FORMAT(26h The data for case number ,I4,3H : ,$) 9000 FORMAT(A1) 9010 FORMAT(A10) 9020 FORMAT(I4,6X,A10) 9030 FORMAT(A30) 9040 FORMAT(A60) RETURN END CC ------------------------------------------------------------------------ CC *JNTGR*: SUBROUTINE FOR TRANSFORMING A CHARACTER INTO AN INTEGER CC ------------------------------------------------------------------------ FUNCTION JNTGR(KAR) CHARACTER*1 KAR IF (KAR.EQ.'0') JNTGR=0 IF (KAR.EQ.'1') JNTGR=1 IF (KAR.EQ.'2') JNTGR=2 IF (KAR.EQ.'3') JNTGR=3 RETURN END CC ------------------------------------------------------------------------ CC *SMISSING*: SUBROUTINE FOR HANDLING MISSING VALUES CC ------------------------------------------------------------------------ SUBROUTINE SMISSING(NVAR,NCAS,MAXP1,MAXN,JCST,X,Y,AW,JNDEX, 1 MVAL,JMISS,JSUBS,VALMS,NSTOP,LAB,JPRT,LUB) DIMENSION X(MAXP1,MAXN),Y(MAXN),AW(MAXN) DIMENSION VALMS(MAXP1) INTEGER JNDEX(MAXN),JSUBS(MAXP1),JMISS(MAXP1) CHARACTER*10 LAB(MAXP1) NVAD=NVAR+1 JHALT=0 JHLT=0 MAXM=(NCAS*4)/5 DO 10 J=1,NVAD JSUBS(J)=0 IF (.NOT.((J.EQ.NVAR.AND.JCST.EQ.1).OR.(JMISS(J).EQ.0)))THEN DO 20 JNC=1,NCAS IF (X(J,JNC).EQ.VALMS(J)) THEN JSUBS(J)=JSUBS(J)+1 JNDEX(JNC)=0 ENDIF 20 CONTINUE IF (JSUBS(J).GT.MAXM) JHALT=JHALT+1 IF (MVAL.EQ.1.AND.JPRT.NE.0) WRITE(LUB,8000) LAB(J),JSUBS(J) ENDIF 10 CONTINUE IF(JHALT.EQ.0) GOTO 30 WRITE(LUB,8010) JHALT,MAXM NSTOP=1 RETURN 30 NCASM=NCAS MAXM=NCAS-MAXM JL=0 JHALT=0 DO 50 JNC=1,NCASM IF (JNDEX(JNC).NE.0) GOTO 65 IF (JHALT.EQ.1.OR.MVAL.EQ.2) GOTO 60 IF (JPRT.NE.0) WRITE(LUB,8020) NVAD JHALT=1 60 JTVAR=0 DO 70 J=1,NVAD IF(((JCST.EQ.1.AND.J.EQ.NVAR).OR.(JMISS(J).EQ.0)).OR. 1 (X(J,JNC).NE.VALMS(J))) GOTO 70 JTVAR=JTVAR+1 JSUBS(JTVAR)=J 70 CONTINUE IF (JTVAR.EQ.0) GOTO 65 IF (MVAL.EQ.1.OR.JHLT.EQ.1) GOTO 80 IF (JPRT.NE.0) WRITE(LUB,8030) JHLT=1 80 IF (MVAL.EQ.1.AND.JPRT.NE.0) THEN IF (JTVAR.GT.10) GOTO 90 WRITE(LUB,8040) JNC,(JSUBS(J),J=1,JTVAR) GOTO 100 90 WRITE(LUB,8040) JNC,(JSUBS(J),J=1,10) WRITE(LUB,8050) (JSUBS(J),J=11,JTVAR) ENDIF 100 IF (MVAL.EQ.2.AND.X(NVAD,JNC).EQ.VALMS(NVAD)) THEN IF (JPRT.NE.0) WRITE(LUB,8060) JNC JTVAR=JTVAR-1 ENDIF IF (JCST.EQ.1) JTVAR=JTVAR+1 IF(MVAL.EQ.2.AND.JTVAR.EQ.NVAR.AND.JPRT.NE.0) WRITE(LUB,8070) JNC IF (.NOT.(MVAL.EQ.2.AND.JTVAR.LT.NVAR.AND.X(NVAD,JNC). 1 NE.VALMS(NVAD))) GOTO 50 65 JL=JL+1 DO 110 J=1,NVAD 110 X(J,JL)=X(J,JNC) JNDEX(JL)=JNC 50 CONTINUE NCAS=JL IF (MVAL.EQ.1) WRITE(LUB,8080) NCAS IF (NCAS.LE.(NVAR*1.75)) THEN WRITE (*,8090) NCAS,NVAR IF (JCST.EQ.1) WRITE(*,8100) NSTOP=1 RETURN ENDIF IF (NCAS.LE.MAXM) THEN WRITE(LUB,8110) NSTOP=1 RETURN ENDIF IF (MVAL.NE.1) THEN DO 150 J=1,NVAR IF ((JCST.EQ.1.AND.J.EQ.NVAR).OR. 1 (JMISS(J).EQ.0)) GOTO 150 JJ=0 NCASM=NCAS JPLUS=0 DO 160 JNC=1,NCAS JPLUS=JPLUS+1 AW(JPLUS)=X(J,JNC) IF (X(J,JNC).EQ.VALMS(J)) THEN NCASM=NCASM-1 JJ=JJ+1 JPLUS=JPLUS-1 Y(JJ)=JNC ENDIF 160 CONTINUE IF (JJ.EQ.0) THEN IF (JPRT.NE.0) WRITE(LUB,8120) LAB(J) ELSE NCASM=NCAS-JJ AMED=AMDAN(AW,MAXN,AW,NCASM) DO 180 J2=1,JJ JY=Y(J2)+0.2 180 X(J,JY)=AMED IF (JPRT.NE.0) WRITE(LUB,8130) LAB(J),JJ DO 190 J2=1,JJ J2Y=Y(J2)+0.2 190 IF (JPRT.NE.0) WRITE(LUB,8140) JNDEX(J2Y) ENDIF 150 CONTINUE ENDIF IF (MVAL.EQ.1) WRITE(*,8080) NCAS 8000 FORMAT(10h Variable ,A10,24h has a missing value for, 1 I4,7h cases.) 8010 FORMAT(15h There are(is) ,I4,29h variable(s),which contain(s)/ 1 24h more than 80 percent (=,I4,28h ) cases with missing value.) 8020 FORMAT(/39h Case has a missing value for variables, 1 18h (variable number ,I4,17h is the response)/1X, 1 4(1H-),25X,9(1H-)) 8030 FORMAT(39h The following cases have been deleted./) 8040 FORMAT(1X,I4,23X,10I4) 8050 FORMAT(28X,10I4) 8060 FORMAT(6h Case ,I4,24h has a missing value for, 1 23h the response variable.) 8070 FORMAT(6h Case ,I4,32h has a missing value for all the, 1 23h explanatory variables.) 8080 FORMAT(/' After treatment of missing values, ',I4, 1 ' cases remain.'/) 8090 FORMAT(/46h Too many coefficients according to the number, 1 10h of cases.//16h Number of cases,8X,2h= ,I5/ 1 26h Number of coefficients = ,I5/25h The number of cases must, 1 37h be twice the number of coefficients.) 8100 FORMAT(31h (Including the constant term!)) 8110 FORMAT(52h More than 80 percent of the cases had to be deleted/ 1 61h because of the missing values. The analysis will be stopped.) 8120 FORMAT(9h Variable ,A10,28h contains no missing values.) 8130 FORMAT(15h The values of ,A10,17h will be replaced, 1 18h by the median for,I5,13h cases namely) 8140 FORMAT(30X,11hcase number,I6) RETURN END CC ------------------------------------------------------------------------ CC SUBROUTINE FOR CALCULATING THE NUMBER OF REPLICATIONS IN THE CC LMS/LTS ALGORITHM CC ------------------------------------------------------------------------ SUBROUTINE SUBREP(NCAS,NVAR,MAXP1,ALGO,NREP,JREPTAB,JREPLOW, 1 JREPHI) INTEGER JREPTAB(MAXP1) INTEGER JREPLOW(MAXP1) INTEGER JREPHI(MAXP1) CHARACTER ALGO,YNOK WRITE(*,8000) NVAR IF (NCAS.LE.JREPLOW(NVAR)) THEN CC-----TAKE ALL COMBINATIONS ALGO='A' NREP=NCOMB(NVAR,NCAS) WRITE(*,8010) ELSE IF (NCAS.GT.JREPHI(NVAR)) THEN CC-----TAKE RANDOM SUBSETS ALGO='R' WRITE(*,8020) NVAR,NCAS,NVAR ELSE CC-----OTHERWISE CHOOSE BETWEEN RANDOM SUBSETS OR ALL COMBINATIONS NREP=NCOMB(NVAR,NCAS) WRITE(*,8030) NREP,NVAR,NCAS IF (NVAR.LE.2) THEN IF (NREP.GE.50000) WRITE(*,8040) ELSE IF (NREP.GE.10000) WRITE(*,8040) ENDIF WRITE(*,8050) JREPTAB(NVAR),NVAR,NREP,NVAR 10 READ(*,9000) YNOK IF (.NOT.(YNOK.EQ.'1'.OR.YNOK.EQ.'2')) THEN WRITE(*,8060) GOTO 10 ENDIF MMM=JNTGR(YNOK) ALGO='R' IF (MMM.EQ.2) ALGO='A' ENDIF IF (ALGO.EQ.'R') THEN CC-----CHOOSE NUMBER OF RANDOM SUBSETS NREP=JREPTAB(NVAR) WRITE(*,8070) NREP 20 READ(*,9000) YNOK IF (YNOK.EQ.'y') YNOK='Y' IF (YNOK.EQ.'n') YNOK='N' IF (.NOT.(YNOK.EQ.'Y'.OR.YNOK.EQ.'N')) THEN WRITE(*,8060) GOTO 20 ENDIF IF (YNOK.EQ.'N') THEN WRITE(*,8080) READ(*,*) NREP ENDIF ENDIF ENDIF 8000 FORMAT(//' The program will use subsets containing', 1 I4,' cases.') 8010 FORMAT(/' All subsets will be considered.') 8020 FORMAT(/' There are over a million subsets', 1 ' of',I4,' cases out of',I5,' cases.'/' Therefore, the', 1 ' program will draw random subsets of',I4,' cases.') 8030 FORMAT(/' There are',I8,' subsets of',I4,' cases', 1 ' out of',I6,' cases.') 8040 FORMAT(/' It would take a lot of time for the', 1 ' program to consider'/' all these subsets !') 8050 FORMAT(/' Do you agree to have the program consider ',I6, 1 /' randomly drawn subsets of',I4,' cases ? (Then', 1 ' type 1)'/' Or do you prefer instead to have the', 1 ' program construct'/' all',I8,' subsets of',I4, 1 ' cases ? (Then type 2) : ',$) 8060 FORMAT(/' Not allowed ! Enter your choice again : ',$) 8070 FORMAT(/' Do you agree with the default choice of', 1 1X,I6,' subsets ? Yes or no : ',$) 8080 FORMAT(/' Please enter the number of subsets you want : ',$) 9000 FORMAT(A1) RETURN END CC --------------------------------------------------------------------- CC *NCOMB* COMPUTES THE NUMBER OF COMBINATIONS OF K POINTS OUT OF N CC --------------------------------------------------------------------- FUNCTION NCOMB(K,N) COMB=1.0 DO 10 J=1,K FACT=(N-J+1.0)/(K-J+1.0) COMB=COMB*FACT 10 CONTINUE NCOMB=INT(COMB+0.5) RETURN END CC ------------------------------------------------------------------ CC *STATIS* : STANDARDIZATION OF THE OBSERVATIONS AND CC CALCULATION OF SOME STATISTICS CC ------------------------------------------------------------------ SUBROUTINE STATIS(X,Y,XMED,XMAD,AW2,WEIGHTS,AW, 1 JNDEX,C,H,JSUBS,NSTOP,NVAR,NVAD,NFAC,NCAS,JCST,JPRT, 1 MAXP1,MAXN,LUB,MAXP,PREC,LAB) DIMENSION X(MAXP1,MAXN),Y(MAXN),XMED(MAXP1),XMAD(MAXP1) DIMENSION AW2(MAXN),WEIGHTS(MAXN),AW(MAXN) DIMENSION C(MAXP,MAXP1) DOUBLE PRECISION H(MAXP,MAXP1) INTEGER JNDEX(MAXN),JSUBS(MAXP1) CHARACTER*10 LAB(MAXP1) AL=NCAS IF (JCST.NE.0) GOTO 60 CC---REGRESSION WITHOUT INTERCEPT DO 50 J=1,NVAD XMED(J)=0.0 DO 10 JNC=1,NCAS 10 AW2(JNC)=ABS(X(J,JNC)) XMAD(J)=AMDAN(AW,MAXN,AW2,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)+AW2(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 CC---REGRESSION WITH INTERCEPT 60 XMED(NVAR)=0.0 XMAD(NVAR)=1.0 DO 120 J=1,NVAD IF(J.EQ.NVAR) GOTO 120 DO 70 JNC=1,NCAS 70 AW2(JNC)=X(J,JNC) XMED(J)=AMDAN(AW,MAXN,AW2,NCAS) DO 80 JNC=1,NCAS 80 AW2(JNC)=ABS(AW2(JNC)-XMED(J)) XMAD(J)=AMDAN(AW,MAXN,AW2,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)+AW2(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 CC---PRINT MEDIANS WRITE(LUB,8010) WRITE(LUB,8020) (LAB(J),J=1,NFAC),LAB(NVAD) IF (NVAD.GT.7) GOTO 130 IF (JCST.EQ.1) WRITE(LUB,8030)(XMED(J),J=1,NFAC),XMED(NVAD) GOTO 140 130 WRITE(LUB,8030)(XMED(J),J=1,6) IF (NFAC.GE.7) WRITE(LUB,8040)(XMED(J),J=7,NFAC),XMED(NVAD) IF (NFAC.LT.7) WRITE(LUB,8040) XMED(NVAD) CC---PRINT DISPERSIONS 140 WRITE(LUB,8050) 150 IF (JCST.EQ.1) WRITE(LUB,8020) (LAB(J),J=1,NFAC),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 160 IF (JCST.EQ.0) WRITE(LUB,8030)(XMAD(J),J=1,NVAD) IF (JCST.EQ.1) WRITE(LUB,8030)(XMAD(J),J=1,NFAC),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.NFAC.GE.7) 1 WRITE(LUB,8040)(XMAD(J),J=7,NFAC),XMAD(NVAD) IF (JCST.EQ.1.AND.NFAC.LT.7) WRITE(LUB,8040) XMAD(NVAD) 170 IF (JPRT.LT.2) GOTO 200 WRITE(LUB,8060) IF (JCST.EQ.1) WRITE(LUB,8025) (LAB(J),J=1,NFAC),LAB(NVAD) IF (JCST.EQ.0) WRITE(LUB,8025) (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) JNDEX(JNC),(X(J,JNC),J=1,NFAC), 1 X(NVAD,JNC) IF (JCST.EQ.0) WRITE(LUB,8070) JNDEX(JNC),(X(J,JNC),J=1,NVAD) GOTO 190 180 WRITE(LUB,8070) JNDEX(JNC),(X(J,JNC),J=1,6) IF (JCST.EQ.1.AND.NFAC.GE.7) 1 WRITE(LUB,8040) (X(J,JNC),J=7,NFAC),X(NVAD,JNC) IF (JCST.EQ.1.AND.NFAC.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,NFAC,JCST,AW,Y,AW2,MAXP1,MAXN, 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,NFAC 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,NFAC),DIAG IF (JCST.EQ.0) WRITE(LUB,8090) LAB(NVAD),(H(JSPA,JSPB), 1 JSPB=1,NVAR),DIAG DO 240 J=1,NVAR 240 JSUBS(J)=J WRITE(LUB,8100) LAB(NVAD) WRITE(LUB,8090) LAB(1),DIAG DO 250 JSPA=1,NFAC 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,NFAC),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(/21h Variable with label ,A10,17h is constant over/ 1 47h all observations hence the program cannot run., 1 /50h Please repeat the analysis without this variable./) 8010 FORMAT(//12h Medians = /) 8015 FORMAT(//33h Dispersion of absolute values = /) 8020 FORMAT(7X,6(A10,1X)/15X,5(A10,1X)) 8025 FORMAT(' i ',6(A10,1X)/15X,5(A10,1X)) 8030 FORMAT(7X,6(F10.4,1X)) 8040 FORMAT(15X,5(F10.4,1X)) 8050 FORMAT(//15h Dispersions = /) 8060 FORMAT(//35h The standardized observations are:/) 8070 FORMAT(I5,2X,6(F10.4,1X)) 8080 FORMAT(//33h Pearson correlation coefficients, 1 22h between the variables/' (the response variable is ', 1 A10,')'/) 8090 FORMAT(1X,A10,(2X,11F6.2)/) 8100 FORMAT(//39h Spearman rank correlation coefficients, 1 22h between the variables/' (the response variable is ', 1 A10,')'/) RETURN END CC ----------------------------------------------------------------- CC *CORR* : CALCULATES THE SPEARMAN RANK CORRELATION COEFFICIENT CC BETWEEN VARIABLES JSPA AND JSPB. CC ----------------------------------------------------------------- SUBROUTINE CORR(X,JSPB,JSPA,NCAS,NFAC,JCST,AW,Y,AW2,MAXP1,MAXN, 1 SPEAR,PEARC,PREC) DIMENSION X(MAXP1,MAXN),AW(MAXN),Y(MAXN),AW2(MAXN) 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.NFAC)) 1 GOTO 5 GOTO 105 5 NTWE=0 DO 10 NJ=1,NCAS AW(NJ)=X(JSPA,NJ) AW2(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=AW2(NJA) AW2(NJA)=AW2(NJ) AW2(NJ)=W 20 CONTINUE DO 40 NJ=1,NCAS NDX=AW2(NJ) IF(NTWE.EQ.0) Y(NDX)=NJ IF(NTWE.NE.0) AW(NDX)=NJ 40 CONTINUE NJ=1 55 TA=1 NJA=AW2(NJ) IF(NTWE.EQ.0) TB=Y(NJA) IF(NTWE.NE.0) TB=AW(NJA) 60 NJA=AW2(NJ) NJAP=AW2(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=AW2(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) AW2(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'. CC ----------------------------------------------------------------- FUNCTION PULL(AW,MAXN,AA,NCAS,K) DIMENSION AA(NCAS),AW(MAXN) CC-----AW IS THE WORKING VECTOR, AA DOES NOT CHANGE 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'. CC ----------------------------------------------------------------- FUNCTION AMDAN(AW,MAXN,AA,NCAS) DIMENSION AA(NCAS),AW(MAXN) CC-----AW IS THE WORKING VECTOR, AA DOES NOT CHANGE AL=NCAS JNDL=INT(AL/2.0) IF (MOD(NCAS,2).EQ.0) THEN AMDAN= (PULL(AW,MAXN,AA,NCAS,JNDL)+ 1 PULL(AW,MAXN,AA,NCAS,JNDL+1))/2.0 ELSE AMDAN= PULL(AW,MAXN,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,MAXN,SSE, 1 FVALUE,PREC,XMAD,XMED,MAXP1,WEIGHTS,NNUL) DIMENSION Y(MAXN),XMAD(MAXP1),XMED(MAXP1),WEIGHTS(MAXN) 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. CC ----------------------------------------------------------------- SUBROUTINE RESTT(GESCF,JABS,JTR,LTSTE,NCAS,NVAD,NZWE, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,MAXP1,MAXN,NOPT,SCAL,QUAN,PREC) DIMENSION GESCF(LTSTE),X(MAXP1,MAXN),Y(MAXN),WEIGHTS(MAXN) DIMENSION RESDU(MAXN),XMED(MAXP1),XMAD(MAXP1) 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/LTS 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 IF (NOPT.EQ.0) RETURN SCAL=SQRT(SCAL/(AVLQS-ANVAR)) 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 40 CONTINUE RETURN END CC ----------------------------------------------------------------- CC *RANGS* : SORTS A VECTOR OF LENGTH 'NCAS'. CC ----------------------------------------------------------------- SUBROUTINE RANGS(AM,LAME) DIMENSION AM(LAME) INTEGER JLV(30),JRV(30) JSS=1 JLV(1)=1 JRV(1)=LAME 10 JNDL=JLV(JSS) JR=JRV(JSS) JSS=JSS-1 20 JNC=JNDL J=JR JTWE=(JNDL+JR)/2 XX=AM(JTWE) 30 IF(AM(JNC).GE.XX)GOTO 40 JNC=JNC+1 GOTO 30 40 IF(XX.GE.AM(J))GOTO 50 J=J-1 GOTO 40 50 IF(JNC.GT.J)GOTO 60 AMM=AM(JNC) AM(JNC)=AM(J) AM(J)=AMM JNC=JNC+1 J=J-1 60 IF (JNC.LE.J) GOTO 30 IF ((J-JNDL).LT.(JR-JNC)) GOTO 80 IF (JNDL.GE.J) GOTO 70 JSS=JSS+1 JLV(JSS)=JNDL JRV(JSS)=J 70 JNDL=JNC GOTO 100 80 IF (JNC.GE.JR) GOTO 90 JSS=JSS+1 JLV(JSS)=JNC JRV(JSS)=JR 90 JR=J 100 IF (JNDL.LT.JR) GOTO 20 IF (JSS.NE.0) GOTO 10 RETURN END CC ------------------------------------------------------------------- CC *RDUAL* : CALCULATES THE RESIDUALS OF ALL CASES CC ------------------------------------------------------------------- SUBROUTINE RDUAL(AA,JKEUS,JDA,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC, 1 JREG,X,Y,RESDU,WEIGHTS,XMED,XMAD,JPLT,AW,JNDEX, 1 MAXP1,MAXN,JHEAD,LAB,MROB) DIMENSION AA(JDA),X(MAXP1,MAXN),Y(MAXN),RESDU(MAXN) DIMENSION WEIGHTS(MAXN) DIMENSION XMED(MAXP1),XMAD(MAXP1),AW(MAXN) INTEGER JNDEX(MAXN) CHARACTER*10 LAB(MAXP1) CHARACTER*60 JHEAD JPL2=1 IF(JPRT.EQ.0) GOTO 10 IF (NVAR.EQ.1.AND.JCST.EQ.1) THEN IF (JKEUS.NE.2) WRITE(LUB,8000) LAB(NVAD) IF (JKEUS.EQ.2) WRITE(LUB,8010) LAB(NVAD) ELSE IF (JKEUS.NE.2) WRITE(LUB,8020) LAB(NVAD),LAB(NVAD) IF (JKEUS.EQ.2) WRITE(LUB,8030) LAB(NVAD),LAB(NVAD) ENDIF 10 DO 20 JNC=1,NCAS IF (.NOT.(NVAR.EQ.1.AND.JCST.EQ.1)) GOTO 40 BBB=Y(JNC)-AA(1) GOTO 60 40 AW(JNC)=0.0 DO 30 J=1,NVAR 30 AW(JNC)=AW(JNC)+AA(J)*(X(J,JNC)*XMAD(J)+XMED(J)) YJ=Y(JNC)*XMAD(NVAD)+XMED(NVAD) BBB=YJ-AW(JNC) 60 IF (AA(NVAD).GT.PREC) GOTO 70 IF (ABS(BBB).LT.PREC) BBB=0.0 RESDU(JNC)=BBB GOTO 50 70 RESDU(JNC)=BBB/AA(NVAD) 50 IF (JPRT.EQ.0) GOTO 20 IF (NVAR.EQ.1.AND.JCST.EQ.1) GOTO 80 IF (AA(NVAD).GT.PREC.AND.JKEUS.NE.2) 1 WRITE(LUB,8040) YJ,AW(JNC),BBB,JNDEX(JNC),RESDU(JNC) IF (AA(NVAD).LE.PREC) 1 WRITE(LUB,8050) YJ,AW(JNC),BBB,JNDEX(JNC) IF (AA(NVAD).GT.PREC.AND.JKEUS.EQ.2) 1 WRITE(LUB,8040) YJ,AW(JNC),BBB,JNDEX(JNC),RESDU(JNC), 1 WEIGHTS(JNC) GOTO 20 80 IF (JKEUS.NE.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC) IF (JKEUS.EQ.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC), 1 WEIGHTS(JNC) 20 CONTINUE IF(JPLT.LT.1) GOTO 90 IF (AA(NVAD).LE.PREC) JPL2=0 CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,JNDEX,MAXN,LUB,JREG, 1 JHEAD,MAXP1,NVAD,LAB,MROB) IF (JPLT.LT.3) GOTO 90 JPLT=2 CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,JNDEX,MAXN,LUB,JREG, 1 JHEAD,MAXP1,NVAD,LAB,MROB) JPLT=3 90 CONTINUE 8000 FORMAT(//9X,8hobserved,12X,'residual',5h i,8h res/s/7X,A10/) 8010 FORMAT(//9X,8hobserved,12X,'residual',5h i,8h res/s,' w'/ 1 7X,A10,33X,5h i/) 8020 FORMAT(//9X,8hobserved,11X,9hestimated,12X,8hresidual,3X,2h i, 1 8h res/s/,7X,A10,10X,A10/) 8030 FORMAT(//9X,8hobserved,11X,9hestimated,12X,'residual',3X,2h i, 1 8h res/s,' w'/7X,A10,10X,A10,33X,5h i/) 8040 FORMAT(1X,F16.5,2(1X,F19.5),I5,1X,F7.2,F6.1) 8050 FORMAT(1X,F16.5,2(1X,F19.5),I5,8H *****) 8060 FORMAT(1X,F16.5,1X,F19.5,I5,1X,F7.2,F6.1) RETURN END CC ---------------------------------------------------------------- CC *LMSLOC* : CALCULATES THE LMS IN THE ONE-DIMENSIONAL CASE. CC W CONTAINS THE ORDERED OBSERVATIONS CC ---------------------------------------------------------------- SUBROUTINE LMSLOC(W,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) DIMENSION W(NCAS),AW2(NCAS) K=JQU-1 VERSC=W(K+1)-W(1) SLUTN=(W(1)+W(K+1))/2.0 AW2(1)=SLUTN BSTD=(W(K+1)-W(1))/2.0 IF (NCAS.EQ.2) RETURN JTEL=1 LNK=NCAS-K DO 10 JNC=2,LNK JPK=JNC+K ABSW=W(JPK)-W(JNC) IF ( ABS(ABSW-VERSC).GT.PREC) GOTO 20 JTEL=JTEL+1 AW2(JTEL)=(W(JNC)+W(JPK))/2.0 GOTO 10 20 IF (ABSW.LT.VERSC) THEN VERSC=ABSW JTEL=1 BSTD=(W(JPK)-W(JNC))/2.0 AW2(1)=(W(JNC)+W(JPK))/2.0 ENDIF 10 CONTINUE CC-----IN CASE OF SEVERAL MINIMA, TAKE LOW MEDIAN SLUTN=AW2(INT(JTEL/2.0+0.5)) BSTD=FACTOR*BSTD RETURN END CC ---------------------------------------------------------------- CC *LTSLOC* : CALCULATES THE LTS IN THE ONE-DIMENSIONAL CASE. CC W CONTAINS THE ORDERED OBSERVATIONS CC ---------------------------------------------------------------- SUBROUTINE LTSLOC(W,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) DIMENSION W(NCAS),AW2(NCAS) SQ=0.0 SQMIN=0.0 JTEL=1 DO 20 JINT=1,NCAS-JQU+1 AW2(JINT)=0.0 DO 10 J=1,JQU AW2(JINT)=AW2(JINT)+W(J+JINT-1) IF (JINT.EQ.1) SQ=SQ+W(J)*W(J) 10 CONTINUE SS1=AW2(JINT)*AW2(JINT)/JQU IF (JINT.EQ.1) THEN SQ=SQ-SS1 SQMIN=SQ ELSE SQ=SQ-W(JINT-1)*W(JINT-1)+ W(JINT+JQU-1)*W(JINT+JQU-1) * -SS1 + SS2 IF (ABS(SQ-SQMIN).GT.PREC) THEN IF(SQ.LT.SQMIN) THEN JTEL=1 SQMIN=SQ AW2(1)=AW2(JINT) ENDIF ELSE JTEL=JTEL+1 AW2(JTEL)=AW2(JINT) ENDIF ENDIF SS2=SS1 20 CONTINUE SLUTN=AW2(INT((JTEL+1)/2))/JQU BSTD=FACTOR*SQRT(SQMIN/JQU) RETURN END CC ----------------------------------------------------------------- CC *TRC* : RETRANSFORMATION OF THE VARIANCE-COVARIANCE MATRIX CC ----------------------------------------------------------------- SUBROUTINE TRC(H,DA,MAXP,MAXP1,NVAR,JCST,NFAC,NVAD, 1 XMED,XMAD) DIMENSION DA(MAXP) DIMENSION XMED(MAXP1),XMAD(MAXP1) DOUBLE PRECISION H(MAXP,MAXP1),XMP2,HNN XMP2=DBLE(XMAD(NVAD))*DBLE(XMAD(NVAD)) IF (JCST.EQ.0) THEN DO 10 J=1,NVAR DO 20 K=1,J H(J,K)=H(J,K)*(XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K)))) 20 CONTINUE DA(J)=DSQRT(H(J,J)) 10 CONTINUE ELSE DO 25 J=1,NVAR H(J,NVAD)=H(J,J) 25 CONTINUE DO 30 J=1,NVAR DO 40 K=1,J H(J,K)=H(J,K)*XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K))) 40 CONTINUE DA(J)=DSQRT(H(J,J)) 30 CONTINUE DO 50 K=1,NFAC H(NVAR,K)=H(K,NVAR)*XMP2/DBLE(XMAD(K)) DO 60 K2=1,NVAR IF (K.EQ.K2) THEN H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K))*XMP2/ 1 (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,NVAD) GOTO 60 ENDIF IF (K.LT.K2) THEN H(NVAR,K)=H(NVAR,K)-(DBLE(XMED(K2))*XMP2)/ 1 (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K,K2) ELSE H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K2))*XMP2/ 1 (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,K) ENDIF 60 CONTINUE 50 CONTINUE H(NVAR,NVAR)=H(NVAR,NVAD)*XMP2 DO 70 K=1,NVAR H(NVAR,NVAR)=H(NVAR,NVAR)+ 1 (DBLE(XMED(K))*DBLE(XMED(K)))*XMP2/ 1 (DBLE(XMAD(K))*DBLE(XMAD(K)))*H(K,NVAD) 70 CONTINUE DO 80 K=1,NVAR IF (K.NE.NVAR) THEN H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/ 1 (DBLE(XMAD(K)))*H(K,NVAR) ELSE H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/ 1 (DBLE(XMAD(K)))*H(NVAR,NVAD) ENDIF 80 CONTINUE DO 90 J=1,NFAC JU=J+1 DO 90 K=JU,NVAR HNN=2.0D0*DBLE(XMED(J))*DBLE(XMED(K))*XMP2 H(NVAR,NVAR)=H(NVAR,NVAR)+HNN/ 1 (DBLE(XMAD(J))*DBLE(XMAD(K)))*H(J,K) 90 CONTINUE DA(NVAR)=DSQRT(H(NVAR,NVAR)) ENDIF RETURN END CC ----------------------------------------------------------------- CC *SCHCV* : PRINTS THE VARIANCE-COVARIANCE MATRIX OF THE CC REGRESSION COEFFICIENTS. CC ----------------------------------------------------------------- SUBROUTINE SCHCV(NVAR,MAXP1,MAXP,H,LUB) DOUBLE PRECISION H(MAXP,MAXP1) WRITE(LUB,9000) DO 10 J=1,NVAR WRITE(LUB,9010) (H(J,K),K=1,J) 10 CONTINUE 9000 FORMAT(/' Variance-covariance matrix of the estimated', 1 ' coefficients :'/) 9010 FORMAT(5(5X,D10.4)) RETURN END CC ---------------------------------------------------------------- CC *RTRAN* : RETRANSFORMATION OF THE REGRESSION COEFFICIENTS. CC ---------------------------------------------------------------- SUBROUTINE RTRAN(NVAR,JCST,NFAC,NVAD,MAXP1,XMED,XMAD,AA,JAL,FCKW) DIMENSION AA(JAL),XMED(MAXP1),XMAD(MAXP1) IF (NVAR.LE.1) THEN AA(1)=AA(1)*XMAD(NVAD)/XMAD(1) GOTO 30 ENDIF DO 10 J=1,NFAC 10 AA(J)=AA(J)*XMAD(NVAD)/XMAD(J) IF (JCST.EQ.0) THEN AA(NVAR)=AA(NVAR)*XMAD(NVAD)/XMAD(NVAR) ELSE AA(NVAR)=AA(NVAR)*XMAD(NVAD) DO 20 J=1,NFAC 20 AA(NVAR)=AA(NVAR)-AA(J)*XMED(J) AA(NVAR)=AA(NVAR)+XMED(NVAD) ENDIF 30 FCKW=FCKW*(XMAD(NVAD)*XMAD(NVAD)) RETURN END CC ----------------------------------------------------------------- CC *RANPN* : RANDOM NUMBER GENERATOR . CC ----------------------------------------------------------------- SUBROUTINE RANPN(NCAS,NVAR,JSUBS,MAXP1,JEY,JLP,MAXRP) INTEGER JSUBS(MAXP1) CC WE PROGRAMMED THIS RANDOM NUMBER GENERATOR OURSELVES BECAUSE CC WE WANTED IT TO BE MACHINE INDEPENDENT. IT SHOULD RUN ON MOST CC COMPUTERS BECAUSE THE LARGEST INTEGER USED IS LESS THAN 2**30. CC THE PERIOD IS 2**16=65536, WHICH IS SUFFICIENT FOR THIS CC APPLICATION. AL=NCAS JLP=JLP+1 IF (JLP.LE.MAXRP) GOTO 10 RETURN 10 DO 20 JNC=1,NVAR 30 JEY=JEY*5761+999 JEYQT=JEY/65536 JEY=JEY-JEYQT*65536 RY=JEY RNULE=RY/65536.0 JRAN=INT(RNULE*AL)+1 IF (JNC.EQ.1) GOTO 50 JNCM=JNC-1 DO 40 J=1,JNCM IF (JSUBS(J).EQ.JRAN) GOTO 30 40 CONTINUE 50 JSUBS(JNC)=JRAN 20 CONTINUE RETURN END CC --------------------------------------------------------------- CC *GENPN* : GENERATES ALL COMBINATIONS OF 'NVAR' INTEGERS CC BETWEEN 1 AND 'NCAS'. CC --------------------------------------------------------------- SUBROUTINE GENPN(NCAS,NVAR,JSUBS) INTEGER JSUBS(NVAR) K=NVAR JSUBS(K)=JSUBS(K)+1 10 IF(K.EQ.1.OR.JSUBS(K).LE.(NCAS-(NVAR-K))) GOTO 100 K=K-1 JSUBS(K)=JSUBS(K)+1 DO 50 I=K+1,NVAR JSUBS(I)=JSUBS(I-1)+1 50 CONTINUE GOTO 10 100 RETURN END CC ------------------------------------------------------------------ CC *EQUAT* : SOLVES A SYSTEM OF LINEAR EQUATIONS. CC ------------------------------------------------------------------ SUBROUTINE EQUAT(AM,MAXP,MAXP1,HVEC,MAXPP1,NA,NB,NERR) DIMENSION AM(MAXP,MAXP1) DOUBLE PRECISION HVEC(MAXPP1),TURN,SWAP,DETER JDM=MAXP DETER=1.0D0 N=NA JMAT=N+NB JNK=0 DO 10 J=1,JMAT JNK=(J-1)*MAXP DO 10 NC=1,MAXP JNK=JNK+1 HVEC(JNK)=AM(NC,J) 10 CONTINUE NZNDE=N-1 LCLPL=-JDM DO 120 JHFD=1,N TURN=0.D0 LCLPL=LCLPL+JDM+1 JDEL=LCLPL+N-JHFD DO 40 JNCB=LCLPL,JDEL IF(DABS(HVEC(JNCB))-DABS(TURN)) 40,40,30 30 TURN=HVEC(JNCB) LDEL=JNCB 40 CONTINUE IF(DABS(TURN).LE.1D-8) GOTO 170 50 IF(LDEL-LCLPL) 60,80,60 60 DETER=-DETER LDEL=LDEL-JDM JNCB=LCLPL-JDM DO 70 JNCC=JHFD,JMAT LDEL=LDEL+JDM JNCB=JNCB+JDM SWAP=HVEC(JNCB) HVEC(JNCB)=HVEC(LDEL) 70 HVEC(LDEL)=SWAP 80 DETER=DETER*TURN IF (JHFD.EQ.N) GOTO 120 TURN=1./TURN JNCB=LCLPL+1 DO 90 JNCC=JNCB,JDEL 90 HVEC(JNCC)=HVEC(JNCC)*TURN JNCD=LCLPL JROW=JHFD+1 DO 110 JNCB=JROW,N JNCD=JNCD+1 JNCE=LCLPL JNCF=JNCD DO 100 JNCC=JROW,JMAT JNCE=JNCE+JDM JNCF=JNCF+JDM 100 HVEC(JNCF)=HVEC(JNCF)-HVEC(JNCE)*HVEC(JNCD) 110 CONTINUE 120 CONTINUE NERR=0 NEQA=N+1 JBEGX=NZNDE*JDM+1 DO 150 JNC=NEQA,JMAT JBEGX=JBEGX+JDM JENDX=JBEGX+N JBEGC=N*JDM+1 JENDC=JBEGC+NZNDE DO 140 JNCB=1,NZNDE JENDX=JENDX-1 JBEGC=JBEGC-JDM JENDC=JENDC-JDM-1 HVEC(JENDX)=HVEC(JENDX)/HVEC(JENDC+1) SWAP=HVEC(JENDX) JNCD=JBEGX-1 DO 130 JNCC=JBEGC,JENDC JNCD=JNCD+1 130 HVEC(JNCD)=HVEC(JNCD)-HVEC(JNCC)*SWAP 140 CONTINUE HVEC(JBEGX)=HVEC(JBEGX)/HVEC(1) 150 CONTINUE JNC=-JDM JBEGX=NZNDE*JDM+1 JENDX=JBEGX+NZNDE DO 160 JNCB=NEQA,JMAT JBEGX=JBEGX+JDM JENDX=JENDX+JDM JNC=JNC+JDM JNCD=JNC DO 160 JNCC=JBEGX,JENDX JNCD=JNCD+1 160 HVEC(JNCD)=HVEC(JNCC) GOTO 180 170 NERR=-1 180 JNK=0 DO 190 J=1,JMAT DO 190 NC=1,MAXP JNK=JNK+1 AM(NC,J)=HVEC(JNK) 190 CONTINUE RETURN END CC ------------------------------------------------------------------ CC *LSL* : CALCULATES THE (REWEIGHTED) LS WHEN 'NVAR' EQUALS 1 CC AND 'JCST'=0. CC ------------------------------------------------------------------ SUBROUTINE LSL(NCAS,MAXN,X,Y,RESDU,A,FCKW,H,MAXP,MAXP1) DOUBLE PRECISION H(MAXP,MAXP1) DIMENSION A(MAXP),X(MAXP1,MAXN),Y(MAXN),RESDU(MAXN) AL=NCAS SXY=0.0 SX2=0.0 DO 10 JNC=1,NCAS SXY=SXY+X(1,JNC)*Y(JNC)*RESDU(JNC) SX2=SX2+(X(1,JNC)**2)*RESDU(JNC) 10 CONTINUE A(1)=SXY/SX2 FCKW=0.0 DO 20 JNC=1,NCAS 20 FCKW=FCKW+( (Y(JNC)-X(1,JNC)*A(1))**2 )*RESDU(JNC) H(1,1)=DBLE((FCKW/(AL-1.0D0))/SX2) RETURN END CC ---------------------------------------------------------------- CC *FCN* : PUTS A ROW OF THE MATRIX X IN A VECTOR. CC ---------------------------------------------------------------- SUBROUTINE FCN(K,M,F,X) DIMENSION F(K),X(M) DO 10 J=1,K F(J)=X(J) 10 CONTINUE RETURN END CC ---------------------------------------------------------------- CC *LSREG* : CALCULATES THE LEAST SQUARES REGRESSION ESTIMATES. CC ---------------------------------------------------------------- SUBROUTINE LSREG(MAXP1,MAXN,MAXP,K,N,M,F,X,Y,W,DA, 1 H,FCKW,HVEC,MAXPP1,JMISS) DIMENSION X(MAXP1,MAXN),F(K),Y(N),W(N),DA(K) DOUBLE PRECISION HVEC(MAXPP1),H(MAXP,MAXP1) DOUBLE PRECISION DFCKW,DFACT DOUBLE PRECISION DWJNC,DYJ,DFKA INTEGER JMISS(MAXP1) KPLUS=K+1 DO 10 JNC=1,K DO 20 J=1,KPLUS H(JNC,J) = 0.D0 20 CONTINUE 10 CONTINUE ANUL=0.0 DO 30 JNC=1,N CALL FCN(K,M,F,X(1,JNC)) DWJNC = DBLE(W(JNC)) ANUL=ANUL+W(JNC) DYJ = DBLE(Y(JNC)) DO 40 KA=1,K DFKA = DBLE(F(KA)) H(KA,K+1) = H(KA,K+1) + DWJNC*DFKA*DYJ DO 50 L=1,KA H(KA,L) = H(KA,L) + DWJNC*DFKA*DBLE(F(L)) 50 CONTINUE 40 CONTINUE 30 CONTINUE DO 60 J=1,K DO 70 JNC=1,J H(JNC,J)=H(J,JNC) 70 CONTINUE 60 CONTINUE CALL MATNV(H,MAXP,MAXP1,HVEC,MAXPP1,K,1,JMISS) MM=K+1 FCKW = QLSRG(K,N,M,MAXP1,MAXN,MAXP,F,X,Y,W,H,MM) DO 80 JNC=1,K F(JNC)=H(JNC,K+1) 80 CONTINUE DFCKW=DBLE(FCKW) ANK=ANUL-K DFACT=DBLE(ANK) DFACT=DFCKW/DFACT DO 90 JNC=1,K DO 100 J=1,K H(JNC,J)=H(JNC,J)*DFACT 100 CONTINUE 90 CONTINUE DO 110 JNC=1,K HDA=H(JNC,JNC) DA(JNC)=SQRT(HDA) 110 CONTINUE RETURN END CC ----------------------------------------------------------------- CC *QLSRG* : EVALUATES THE OBJECTIVE FUNCTION FOR LS REGRESSION. CC ----------------------------------------------------------------- FUNCTION QLSRG(K,N,M,MAXP1,MAXN,MAXP,F,X,Y,W,H,MM) DIMENSION F(K),X(MAXP1,MAXN),Y(N),W(N) DOUBLE PRECISION Q,HSUM,H(MAXP,MAXP1) Q=0.D0 DO 30 JNC=1,N CALL FCN(K,M,F,X(1,JNC)) HSUM=0.D0 DO 20 JNCB=1,K 20 HSUM=H(JNCB,MM)*F(JNCB)+HSUM 30 Q=(HSUM-Y(JNC))*(HSUM-Y(JNC))*W(JNC)+Q QLSRG = Q RETURN END CC ----------------------------------------------------------------- CC *MATNV* : PERFORMS A MATRIX INVERSION. CC ----------------------------------------------------------------- SUBROUTINE MATNV(AM,MAXP,MAXP1,HVEC,MAXPP1,NA,NB,JMISS) DOUBLE PRECISION DETER,TURN,SWAP DOUBLE PRECISION HVEC(MAXPP1),AM(MAXP,MAXP1) INTEGER JMISS(MAXP1) DETER=1.0D0 N=NA NPNB=N+NB JNK=0 DO 10 J=1,NPNB JNK=(J-1)*MAXP DO 10 NC=1,MAXP JNK=JNK+1 HVEC(JNK)=AM(NC,J) 10 CONTINUE JDM=MAXP NMA=N-1 JDELC=1-JDM DO 130 JHFD=1,N TURN=0.0D0 JDELC=JDELC+JDM JDLA=JDELC+JHFD-1 JDLB=JDELC +NMA DO 40 JNCB=JDLA,JDLB IF(DABS(HVEC(JNCB))-DABS(TURN)) 40,40,30 30 TURN=HVEC(JNCB) LDEL=JNCB 40 CONTINUE IF(TURN) 50,180,50 50 JPAAL=LDEL-JDELC+1 JMISS(JHFD)=JPAAL IF(JPAAL-JHFD) 80,80,60 60 DETER=-DETER JPAAL=JPAAL-JDM JNCD=JHFD-JDM DO 70 JNC=1,NPNB JPAAL=JPAAL+JDM JNCD=JNCD+JDM SWAP=HVEC(JNCD) HVEC(JNCD)=HVEC(JPAAL) 70 HVEC(JPAAL)=SWAP 80 DETER=DETER*TURN TURN=1./TURN JNCD=JDELC+NMA DO 90 JNC=JDELC,JNCD 90 HVEC(JNC)=-HVEC(JNC)*TURN HVEC(JDLA)=TURN JNCB=JHFD-JDM JPAAL=1-JDM DO 120 JNC=1,NPNB JPAAL=JPAAL+JDM JNCB=JNCB+JDM IF(JNC-JHFD) 100,120,100 100 JCL=JPAAL+NMA SWAP=HVEC(JNCB) JNCD=JDELC-1 DO 110 JNCC=JPAAL,JCL JNCD=JNCD+1 110 HVEC(JNCC)=HVEC(JNCC)+SWAP*HVEC(JNCD) HVEC(JNCB)=SWAP*TURN 120 CONTINUE 130 CONTINUE DO 160 JNCB=1,N JHFD=N+1-JNCB LDEL=JMISS(JHFD) IF(LDEL-JHFD) 140,160,140 140 JPAAL=(LDEL-1)*JDM+1 JCL=JPAAL+NMA JDELC=(JHFD-1)*JDM+1-JPAAL DO 150 JNCC=JPAAL,JCL JNCD=JNCC+JDELC SWAP=HVEC(JNCC) HVEC(JNCC)=HVEC(JNCD) 150 HVEC(JNCD)=SWAP 160 CONTINUE GOTO 180 180 JNK=0 DO 190 J=1,NPNB DO 190 NC=1,MAXP JNK=JNK+1 AM(NC,J)=HVEC(JNK) 190 CONTINUE RETURN END CC ------------------------------------------------------------------- CC *LCAT* : CALCULATES LOCATION ESTIMATES. CC ------------------------------------------------------------------- SUBROUTINE LCAT(NCAS,NVAR,JCST,JPRT,NVAD,X,Y,RESDU,WEIGHTS,PREC, 1 MROB,XMED,XMAD,NZWE,AVW,JPLT,AW,JNDEX,MAXP1,MAXN,MVAL,LUA,LUB, 1 LUC,JREG,JHEAD,FNAMEA,FNAMEB,FNAMEC,YNSAVE,LAB,JFMT,JVARS,YN, 1 JPLACE,AW2,FACLMS,FACLTS) DIMENSION A(11),X(MAXP1,MAXN),Y(MAXN),RESDU(MAXN) DIMENSION XMED(MAXP1),XMAD(MAXP1) DIMENSION WEIGHTS(MAXN),AW(MAXN),AW2(MAXN) DIMENSION FACLMS(11),FACLTS(11) INTEGER JPLACE(MAXP1),JNDEX(MAXN) CHARACTER YN,YNSAVE CHARACTER*10 LAB(MAXP1) CHARACTER*30 FNAMEA,FNAMEB,FNAMEC CHARACTER*60 JFMT,JHEAD CHARACTER*3 NAME NAME=' ' IF (MROB.EQ.1) THEN NAME='LMS' ELSE NAME='LTS' ENDIF WRITE(LUB,8000) WRITE(LUB,8005) JHEAD WRITE(LUB,8010) NCAS IF (MROB.EQ.1) WRITE(LUB,8015) IF (MROB.EQ.2) WRITE(LUB,8016) IF (FNAMEA.NE.'CON') WRITE(LUB,8017) FNAMEA WRITE(LUB,8019) FNAMEB IF (JPRT.EQ.0) WRITE(LUB,8020) IF (JPRT.EQ.1) WRITE(LUB,8040) IF (JPLT.EQ.0) WRITE(LUB,8045) IF (JPLT.NE.0) WRITE(LUB,8050) IF (MVAL.EQ.0) WRITE(LUB,8055) IF (MVAL.EQ.1) WRITE(LUB,8056) DO 5 J=1,11 5 A(J)=0.0 IF (NCAS.LE.2) THEN WRITE(*,8070) NCAS STOP ENDIF IF (FNAMEA.EQ.'CON') WRITE(*,8060) DO 40 JNC=1,NCAS JNDEX(JNC)=JNC IF (FNAMEA.EQ.'CON'.AND.JNC.EQ.1) WRITE(*,8080) JNC IF (FNAMEA.EQ.'CON'.AND.JNC.NE.1) WRITE(*,8090) JNC IF (YN.NE.'Y') GOTO 30 READ(LUA,*) (AW(J),J=1,JVARS) JH=JPLACE(NVAD) Y(JNC)=AW(JH) IF (YNSAVE.EQ.'Y') WRITE(LUC,*) Y(JNC) GOTO 40 30 READ(LUA,JFMT) (AW(J),J=1,JVARS) JH=JPLACE(NVAD) Y(JNC)=AW(JH) IF (YNSAVE.EQ.'Y') WRITE(LUC,JFMT) Y(JNC) 40 X(1,JNC)=Y(JNC) IF (YNSAVE.EQ.'Y') CLOSE(LUC,STATUS='KEEP') IF (MVAL.NE.0) CALL SMISLOC(NCAS,MAXP1,MAXN,X,Y,JNDEX,NSTOP) IF (NCAS.LE.2) THEN WRITE(*,8070) NCAS STOP ENDIF IF (NSTOP.EQ.1) STOP IF (MVAL.NE.0) WRITE(LUB,8095) NCAS AL=NCAS ANVAR=NVAR CALL JQUANT(NCAS,NVAR,MROB,JQU,FACTOR,FACLMS,FACLTS) IF (JPRT.EQ.0) GOTO 60 WRITE(LUB,8100) DO 50 JNC=1,NCAS 50 WRITE(LUB,8110) JNC,Y(JNC) 60 XMED(1)=AMDAN(AW,MAXN,Y,NCAS) DO 70 JNC=1,NCAS RESDU(JNC)=ABS(Y(JNC)-XMED(1)) 70 A(1)=A(1)+Y(JNC) XMAD(1)=AMDAN(AW,MAXN,RESDU,NCAS) XMAD(2)=XMAD(1)*1.4826 WRITE(LUB,8120) XMED(1),XMAD(2) IF (ABS(XMAD(1)).LE.1.0E-12) THEN WRITE(LUB,8125) STOP ENDIF WRITE(LUB,8130) A(1)=A(1)/AL DO 80 JNC=1,NCAS 80 A(2)=A(2)+(Y(JNC)-A(1))**2 A(2)=SQRT(A(2)/(AL-1.0D0)) WRITE(LUB,8140) A(1),A(2) JREG=1 CALL RDUAL(A,0,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,JPLT,AW,JNDEX, 1 MAXP1,MAXN,JHEAD,LAB,MROB) WRITE(LUB,8130) CALL RANGS(Y,NCAS) IF (MROB.EQ.1) THEN CALL LMSLOC(Y,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) FINITF=-10.0/(AL-ANVAR)*(JQU/AL)+(10.0/(AL-ANVAR)+1.0) BSTD=BSTD*FINITF WRITE(LUB,8152) JQU ELSE CALL LTSLOC(Y,NCAS,JQU,SLUTN,BSTD,PREC,AW2,FACTOR) WRITE(LUB,8157) JQU ENDIF JBREAK=NBREAK(JQU,NCAS,NVAR) JDEFAUL=(NCAS+NVAR+1)/2 IF (JBREAK.EQ.NBREAK(JDEFAUL,NCAS,NVAR)) THEN WRITE(LUB,8165) JBREAK,NCAS ELSE WRITE(LUB,8160) JBREAK ENDIF WRITE(LUB,8150) NAME,SLUTN,NAME,BSTD A(1)=SLUTN A(2)=BSTD CC-----DETERMINATION OF THE WEIGHTS NZWE=0 DO 100 JNC=1,NCAS Y(JNC)=X(1,JNC) RESDU(JNC)=(Y(JNC)-A(1))/A(2) IF (ABS(RESDU(JNC)).LT.2.5) THEN WEIGHTS(JNC)=1.0 NZWE=NZWE+1 ELSE WEIGHTS(JNC)=0.0 ENDIF 100 CONTINUE AVW=(NZWE*1.0)/AL AL=NZWE CC-----FINAL SCALE ESTIMATE A(2)=0 DO 105 JNC=1,NCAS 105 A(2)=A(2)+((Y(JNC)-A(1))**2)*WEIGHTS(JNC) A(2)=SQRT(A(2)/(AL-1.0)) WRITE(LUB,8155) NAME,A(2) JREG=2 CALL RDUAL(A,1,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,JPLT,AW,JNDEX, 1 MAXP1,MAXN,JHEAD,LAB,MROB) WRITE(LUB,8135) IF (NZWE.LE.2) THEN WRITE(*,8170) NZWE STOP ENDIF WRITE(LUB,8290) NCAS,JQU IF (MROB.EQ.1) THEN WRITE(LUB,8300) ELSE WRITE(LUB,8305) ENDIF WRITE(LUB,8310) NAME,NAME WRITE(LUB,8130) A(1)=0.0 A(2)=0.0 DO 110 JNC=1,NCAS 110 A(1)=A(1)+Y(JNC)*WEIGHTS(JNC) A(1)=A(1)/AL DO 120 JNC=1,NCAS 120 A(2)=A(2)+((Y(JNC)-A(1))**2)*WEIGHTS(JNC) A(2)=SQRT(A(2)/(AL-1.0)) WRITE(LUB,8180) A(1),A(2) WRITE(LUB,8190) NZWE WRITE(LUB,8200) AVW JREG=3 CALL RDUAL(A,2,NVAD,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,JREG, 1 X,Y,RESDU,WEIGHTS,XMED,XMAD,JPLT,AW,JNDEX, 1 MAXP1,MAXN,JHEAD,LAB,MROB) WRITE (LUB,8130) WRITE(*,8210) IF (YNSAVE.EQ.'Y') WRITE(*,8220) FNAMEC IF (FNAMEA.NE.'CON') WRITE(*,8230) FNAMEA IF (FNAMEB.NE.'CON'.AND.FNAMEB.NE.'PRN') WRITE(*,8240) FNAMEB IF (JPRT.NE.10) STOP 8000 FORMAT(1X, 1 '*******************************************' 1 /' * ROBUST ESTIMATION OF LOCATION AND SCALE *'/ 1 ' *******************************************'///) 8005 FORMAT(/' Title: ',A60) 8010 FORMAT(' Number of cases in the data set: ',I5) 8015 FORMAT(' You have chosen the least median of squares ', 1 '(LMS) method.') 8016 FORMAT(' You have chosen the least trimmed squares ', 1 '(LTS) method.') 8017 FORMAT(' Your data reside in the file: ',A30) 8019 FORMAT(' This output is written in the file: ',A30) 8020 FORMAT(' You have chosen small output.') 8040 FORMAT(' You have chosen large output.') 8045 FORMAT(' You have chosen not to plot the residuals.') 8050 FORMAT(' You have chosen to plot residuals versus their index.') 8055 FORMAT(' The data are assumed not to have missing values.') 8056 FORMAT(' Any case with a missing value will be deleted.') 8070 FORMAT(16h There are only ,I4,25h cases. The analysis must, 1 12h be stopped./) 8060 FORMAT(//24h Please enter your data./) 8080 FORMAT(1X,26h The data for case number ,I4,3H : ,$) 8090 FORMAT(26h The data for case number ,I4,3H : ,$) 8095 FORMAT(/' After treatment of missing values, ',I4, 1 ' cases remain.') 8100 FORMAT(/22h The observations are:/) 8110 FORMAT(I5,2X,F10.4) 8120 FORMAT(//11h The median,8X,2H= ,F11.4,2X, 1 19h The MAD (x 1.4826),10X,2H= ,F11.4/) 8125 FORMAT(/38h More than half of the data are equal.//) 8130 FORMAT(/1X,78(1H*)/) 8135 FORMAT(/1X,78(1H-)/) 8140 FORMAT(//9h The mean,10X,2H= ,F11.4,2X, 1 19h Standard deviation,10X,2H= ,F11.4//) 8150 FORMAT(//1X,A3,' location = ',F11.4,2X, 1 ' Preliminary ',A3,' scale = ',F11.4) 8155 FORMAT(34X,' Final ',A3,' scale estimate = ',F11.4/) 8152 FORMAT(' The program minimizes the ',I4,'th ordered squared', 1 ' residual,') 8157 FORMAT(' The program minimizes the sum of the ',I4,' smallest', 1 ' squared residuals,') 8160 FORMAT(' hence its breakdown value is ',I2,'%.') 8165 FORMAT(' hence its breakdown value is ',I2,'% (= highest', 1 ' possible value '/' for n = ',I4,' cases).') 8290 FORMAT(/' Corresponding formulas: (for this data set ', 1 'n = ',I5,', h = ',I5,')') 8300 FORMAT(// 1 ' LMS objective = obj(y ) = minimum |y - theta|'/ 1 ' i theta i h:n') 8305 FORMAT(// 1 ' | 1 h ', 1 ' 2 |'/ 1 ' LTS objective = obj(y ) = minimum squareroot | - sum', 1 ' |y - theta| |'/ 1 ' i theta | h j=1', 1 ' i j:n |') 8310 FORMAT(// 1 ' Preliminary ',A3,' scale = constant*obj(y )'/ 1 ' i '/// 1 ' | n 2 |'/ 1 ' | sum w *r |'/ 1 ' | i=1 i i |'/ 1 ' Final ',A3,' scale estimate = squareroot | -------------- |'/ 1 ' | n |'/ 1 ' | (sum w ) - 1 |'/ 1 ' | i=1 i |') 8170 FORMAT(/16h There are only ,I4,28h cases with non-zero weight./ 1 30h The analysis must be stopped./) 8180 FORMAT(/' Weighted mean = ',F11.4,2X, 1 ' Weighted standard deviation = ',F11.4/) 8190 FORMAT(/10h There are,2X,I4,30h points with non-zero weight. ) 8200 FORMAT(/' Average weight = ',F8.6/) 8210 FORMAT(////41h The run has been executed successfully. //) 8220 FORMAT(/' The data are saved in the file: ',A30) 8230 FORMAT(' The data has been read from the file: ',A30) 8240 FORMAT(' The output has been written in the file: ',A30) RETURN END CC ---------------------------------------------------------------- CC *SMISLOC* : HANDLING OF MISSING VALUES IN THE CASE OF LOCATION CC ---------------------------------------------------------------- SUBROUTINE SMISLOC(NCAS,MAXP1,MAXN,X,Y,JNDEX,NSTOP) DIMENSION X(MAXP1,MAXN),Y(MAXN) INTEGER JNDEX(MAXN) MAXM=(NCAS*4)/5 JND=0 JCAS=0 CC-----GIVE THE MISSING VALUE CODE WRITE(*,8000) 5 READ(*,*,ERR=5)CODE DO 10 JNC=1,NCAS IF (Y(JNC).EQ.CODE) THEN JND=JND+1 ELSE JCAS=JCAS+1 X(1,JCAS)=Y(JNC) JNDEX(JCAS)=JNC ENDIF 10 CONTINUE NCAS=JCAS WRITE(*,8010) NCAS IF (JND.GT.MAXM) THEN WRITE (*,8020) NSTOP=1 RETURN ENDIF DO 20 JNC=1,NCAS Y(JNC)=X(1,JNC) 20 CONTINUE 8000 FORMAT(/46h Please enter the value of this variable which/ 1 51h has to be interpreted as the missing value code : ,$) 8010 FORMAT(/' After treatment of missing values,', 1 I4,' cases remain.') 8020 FORMAT(/52h More than 80 percent of the cases had to be deleted/ 1 61h because of the missing values. The analysis will be stopped.) RETURN END CC ------------------------------------------------------------------- CC *MOVE* : MOVES A CHARACTER FROM LEFT TO RIGHT CC (USED FOR ADJUSTING VARIABLE LABELS TO THE RIGHT) CC ------------------------------------------------------------------- SUBROUTINE MOVE(KARAK) CHARACTER KARAK(10) NB=0 DO 10 J=1,10 JPL=11-J IF (KARAK(JPL).EQ.' ') NB=NB+1 IF((JPL.NE.1.AND.KARAK(JPL-1).NE.' ').OR.(NB.EQ.10))GOTO 20 10 CONTINUE 20 IF(NB.EQ.10.OR.NB.EQ.0) RETURN NBM=10-NB DO 30 J=1,NBM JPB=NBM+1-J JPL=11-J KARAK(JPL)=KARAK(JPB) KARAK(JPB)=' ' 30 CONTINUE RETURN END CC --------------------------------------------------------------------- CC *PVAL* : THIS IS A FORTRAN VERSION OF THE LACKRITZ ALGORITHM CC (AMERICAN STATISTICIAN 1984, VOL. 38, P.312-314) CC FOR CALCULATING THE ONE-SIDED TAIL AREA OF AN F-DISTRIBUTION CC WITH N1 AND N2 DEGREES OF FREEDOM. IMPLEMENTED BY CC J. VAN SOEST AND B. VAN ZOMEREN, JUNE 1986. CC --------------------------------------------------------------------- FUNCTION PVAL(F,N1,N2) DOUBLE PRECISION F,PI DOUBLE PRECISION PVAL,SUM,T,TINV,T1,TINV1,DSUM,T5,A,B LOGICAL ODD DATA PI/3.14159265/ IF(F.GT.(0.000001))GOTO 5 PVAL=1.0D0 RETURN 5 T=N1*F/N2 TINV=1/T T1=T+1 TINV1=TINV+1 IF (.NOT.ODD(N1) .AND. .NOT. ODD(N2)) THEN M1=N1/2 M2=N2/2 PVAL=SUM(M1-1,DBLE(M2-1),T/T1)/DEXP(DBLE(M2)*DLOG(T1)) ENDIF IF (.NOT.ODD(N1).AND.ODD(N2)) THEN M1=N1/2 M2=(N2-1)/2 A=SUM(M1-1,DBLE(M2-0.5),T/T1) PVAL=A/DEXP(DBLE(M2+0.5)*DLOG(T1)) ENDIF IF (ODD(N1).AND..NOT.ODD(N2)) THEN M1=(N1-1)/2 M2=N2/2 A=SUM(M2-1,DBLE(M1-0.5),TINV/TINV1) PVAL=1-A/DEXP(DBLE(M1+0.5)*DLOG(TINV1)) ENDIF IF (ODD(N1).AND.ODD(N2)) THEN M1=(N1-1)/2 M2=(N2-1)/2 A=0.0 IF (M1.GT.0) THEN DSUM=2*(DSQRT(T))/(PI*T1) DO 10 I=1,M2 DSUM=DSUM*DBLE(I)/(T1*DBLE(I-0.5)) 10 CONTINUE A=DSUM DO 20 I=2,M1 DSUM=DSUM*T*2*DBLE(M2+I-1)/(T1*DBLE(2*I-1)) A=A+DSUM 20 CONTINUE ENDIF B=0.0 IF (M2.GT.0) THEN DSUM=DSQRT(TINV)*2/(PI*TINV1) B=DSUM T5=2*TINV/TINV1 DO 30 I=2,M2 DSUM=DSUM*T5*DBLE(I-1)/DBLE(2*I-1) B=B+DSUM 30 CONTINUE ENDIF PVAL=A-B+(2/PI)*DATAN(DSQRT(TINV)) ENDIF RETURN END CC ------------------------------------------------------------------- CC *PTVAL* : FUNCTION (BASED ON THE FUNCTION PVAL) FOR CALCULATING CC THE TWO-SIDED TAIL AREA OF A STUDENT DISTRIBUTION CC WITH N DEGREES OF FREEDOM CC ------------------------------------------------------------------- FUNCTION PTVAL(TVAL,N) DOUBLE PRECISION PVAL,PTVAL DOUBLE PRECISION TVAL PTVAL=PVAL(TVAL*TVAL,1,N) RETURN END CC ------------------------------------------------------------------- CC *SUM* : AUXILIARY FUNCTION (FOR FUNCTION PVAL) CC ------------------------------------------------------------------- FUNCTION SUM(N,F1,F2) DOUBLE PRECISION F1,F2,DSUM,TSUM,SUM DSUM=1 TSUM=1 IF (N.NE.0) THEN DO 10 I=1,N DSUM=DSUM*(F1+I)*F2/I TSUM=TSUM+DSUM 10 CONTINUE ENDIF SUM=TSUM RETURN END CC ------------------------------------------------------------------- CC *ODD* : AUXILIARY FUNCTION (FOR FUNCTION PVAL) CC ------------------------------------------------------------------- FUNCTION ODD(N) LOGICAL ODD ODD=.TRUE. IF(2*(N/2).EQ.N) ODD=.FALSE. RETURN END CC ------------------------------------------------------------------- CC *GRAF* : PERFORMS RESIDUAL PLOTS AND/OR A PLOT OF THE CC OBSERVATIONS IN THE TWO-DIMENSIONAL CASE. CC ------------------------------------------------------------------- SUBROUTINE GRAF(X,Y,NCAS,JPLXY,JPL2,JPLT,JNDEX,MAXN,LUB,JREG, 1 JHEAD,MAXP1,NVAD,LAB,MROB) DIMENSION X(NCAS),Y(NCAS) INTEGER JNDEX(MAXN),MGRAF(41,51) CHARACTER JDRAW(51),NUM(14) CHARACTER*10 LAB(MAXP1) CHARACTER*60 JHEAD LOGICAL ODD NUM(1)=' ' NUM(2)='1' NUM(3)='2' NUM(4)='3' NUM(5)='4' NUM(6)='5' NUM(7)='6' NUM(8)='7' NUM(9)='8' NUM(10)='9' NUM(11)='*' NUM(12)='-' NUM(13)='+' NUM(14)=' ' JXSCA=51 JYSCA=41 AFR=10.0E-02 DO 10 J=1,JXSCA DO 10 K=1,JYSCA 10 MGRAF(K,J)=0 IF (JPLT.EQ.2.AND.JPLXY.NE.1) X(1)=JNDEX(1) XMIN=X(1) XMAX=X(1) YMAX2=2.5 YKLE2=-2.5 YKLE=Y(1) YMAX=Y(1) DO 20 JNC=1,NCAS IF (JPLT.EQ.2.AND.JPLXY.NE.1) X(JNC)=JNDEX(JNC) IF (X(JNC).GT.XMAX) XMAX=X(JNC) IF (X(JNC).LT.XMIN) XMIN=X(JNC) IF (Y(JNC).GT.YMAX) YMAX=Y(JNC) IF (Y(JNC).LT.YKLE) YKLE=Y(JNC) 20 CONTINUE IF((YMAX-YKLE).LE.AFR.AND.JPLXY.EQ.1)THEN YKLE=YKLE-1.0 YMAX=YMAX+1.0 ENDIF IF (NVAD.EQ.2.AND.JPLXY.EQ.1) THEN IF(0.0.GT.XMAX)XMAX=0.0 IF(0.0.LT.XMIN)XMIN=0.0 IF(0.0.GT.YMAX)YMAX=0.0 IF(0.0.LT.YKLE)YKLE=0.0 ENDIF IF (YMAX.GT.2.5.OR.JPLXY.EQ.1.OR.(JPL2.EQ.0.AND.YKLE.NE.YMAX)) 1 YMAX2=YMAX IF (YKLE.LT.(-2.5).OR.JPLXY.EQ.1.OR.(JPL2.EQ.0.AND.YKLE.NE.YMAX)) 1 YKLE2=YKLE XSCA=JXSCA YSCA=JYSCA IF((XMAX-XMIN).LE.AFR)THEN XMIN=XMIN-1.0 XMAX=XMAX+1.0 ENDIF A=(XSCA-1.0)/(XMAX-XMIN) B=1.0-A*XMIN C=(YSCA-1.0)/(YMAX2-YKLE2) D=1.0-C*YKLE2 IF (JPLXY.EQ.1) GOTO 50 JNUL=INT(YSCA-D) IF ((JNUL.LE.1.OR.ABS(YMAX2).LT.AFR).AND.JPL2.NE.0) THEN J2P=1 JNUL=2 J2M=3 GOTO 50 ENDIF IF ((JNUL.GE.JYSCA.OR.ABS(YKLE2).LT.AFR).AND.JPL2.NE.0) THEN J2P=JYSCA-2 JNUL=JYSCA-1 J2M=JYSCA GOTO 50 ENDIF IF (JPL2.EQ.0) THEN IF (JNUL.LE.1.OR.ABS(YMAX2).LE.AFR) JNUL=1 IF (JNUL.GT.JYSCA.OR.ABS(YKLE2).LE.AFR) JNUL=JYSCA ENDIF IF ((JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 50 J2P=INT(YSCA-(2.5*C+D)) IF (JNUL.EQ.J2P) J2P=JNUL-1 IF (ABS(YMAX2-2.5).LT.AFR) J2P=1 IF (J2P.LE.0) THEN IF (ABS(YMAX2-2.5).LT.AFR) THEN J2P=1 ELSE J2P=2 ENDIF ENDIF IF (J2P.GT.JYSCA) J2P=JYSCA-2 IF (ABS(YKLE2-2.5).LT.AFR) J2P=JYSCA-2 J2M=2*JNUL-J2P IF (J2M.LE.0.OR.ABS(YMAX2+2.5).LT.AFR) THEN J2P=1 JNUL=2 J2M=3 GOTO 50 ENDIF IF (J2M.GT.JYSCA.OR.ABS(YKLE2+2.5).LT.AFR) THEN J2M=JYSCA JS=J2M+J2P IF (ODD(JS)) THEN J2P=J2P+1 ENDIF JNUL=(J2M+J2P)/2 ENDIF 50 DO 110 JNC=1,NCAS XJNC=X(JNC) YJNC=Y(JNC) 60 JX=INT(XJNC*A+B) IF (JX.LE.0) JX=1 IF (JX.GT.JXSCA) JX=JXSCA IF(.NOT.(ABS(YJNC-YMAX2).LT.AFR.OR.ABS(YJNC-YKLE2).LT.AFR)) 1 GOTO 70 IF (ABS(YJNC-YMAX2).LT.AFR) JEY=1 IF (ABS(YJNC-YKLE2).LT.AFR) JEY=JYSCA IF (ABS(YJNC).LT.AFR.AND.JPLXY.EQ.0) JEY=JNUL GOTO 90 70 JEY=INT(YSCA-(YJNC*C+D)) IF (JPLXY.NE.1) THEN IF (ABS(YJNC).LT.AFR) JEY=JNUL IF (JPL2.EQ.0) GOTO 80 IF (JEY.EQ.J2P.AND.YJNC.LT.(2.5-AFR)) JEY=J2P+1 IF (JEY.EQ.J2P.AND.YJNC.GT.(2.5+AFR)) JEY=J2P-1 IF (JEY.EQ.J2M.AND.YJNC.LT.(-2.5-AFR)) JEY=J2M+1 IF (JEY.EQ.J2M.AND.YJNC.GT.(-2.5+AFR)) JEY=J2M-1 IF (ABS(YJNC-2.5).LT.AFR.OR.(JEY.GT.J2P.AND.YJNC.GT.2.5) 1 .OR.(JEY.LT.J2P.AND.YJNC.LT.2.5)) JEY=J2P IF (ABS(YJNC+2.5).LT.AFR.OR.(JEY.GT.J2M.AND.YJNC.GT.-2.5) 1 .OR.(JEY.LT.J2M.AND.YJNC.LT.-2.5)) JEY=J2M ENDIF 80 IF (JEY.LE.0) JEY=1 IF (JEY.GT.JYSCA) JEY=JYSCA 90 MGRAF(JEY,JX)=MGRAF(JEY,JX)+1 IF(NUM(14).NE.'0') GOTO 100 MGRAF(JEY,JX)=-20 GOTO 120 100 IF (.NOT.(JNC.EQ.NCAS.AND.JPLXY.EQ.1.AND.NVAD.EQ.2)) GOTO 110 XJNC=0.0 YJNC=0.0 NUM(14)='0' GOTO 60 110 CONTINUE 120 IF (JREG.EQ.0) WRITE(LUB,8000) JHEAD IF (JREG.EQ.1) WRITE(LUB,8010) JHEAD IF (JREG.EQ.2.AND.MROB.EQ.1) WRITE(LUB,8020) JHEAD IF (JREG.EQ.2.AND.MROB.EQ.2) WRITE(LUB,8025) JHEAD IF (JREG.EQ.3) WRITE(LUB,8030) JHEAD IF (JPLXY.EQ.0.AND.JPL2.NE.0) WRITE(LUB,8050) IF (JPLXY.EQ.0.AND.JPL2.EQ.0) WRITE(LUB,8060) IF (JPLXY.EQ.1) WRITE(LUB,8070) LAB(NVAD) WRITE(LUB,8075) DO 240 J=1,JYSCA IF (JPLXY.EQ.1) GOTO 160 IF(J.NE.JNUL) GOTO 140 DO 130 K=1,JXSCA MJK=MGRAF(J,K)+1 IF (MJK.EQ.1) JDRAW(K)=NUM(12) IF (MJK.GT.10) JDRAW(K)=NUM(11) IF (MJK.GE.2.AND.MJK.LE.10) JDRAW(K)=NUM(MJK) 130 CONTINUE GOTO 180 140 IF ((JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 160 IF (J.NE.J2P.AND.J.NE.J2M) GOTO 160 DO 150 K=1,JXSCA MJK=MGRAF(J,K)+1 IF (MJK.EQ.1) JDRAW(K)=NUM(13) IF (MJK.GT.10) JDRAW(K)=NUM(11) IF (MJK.GE.2.AND.MJK.LE.10) JDRAW(K)=NUM(MJK) 150 CONTINUE GOTO 180 160 DO 170 K=1,JXSCA MJK=MGRAF(J,K)+1 IF(MJK.GT.10) JDRAW(K)=NUM(11) IF(MJK.LE.10) JDRAW(K)=NUM(MJK) IF(MJK.LT.0) JDRAW(K)=NUM(14) 170 CONTINUE 180 JREST=MOD(J,5) IF (JPLXY.EQ.1) GOTO 190 IF ((JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 190 IF (J.NE.J2P) GOTO 190 IF (JREST.EQ.1) WRITE(LUB,8080)(JDRAW(K),K=1,JXSCA) IF (JREST.NE.1) WRITE(LUB,8085)(JDRAW(K),K=1,JXSCA) GOTO 240 190 IF (J.NE.1) GOTO 200 WRITE(LUB,8130) YMAX2,(JDRAW(K),K=1,JXSCA) GOTO 240 200 IF (JPLXY.EQ.1) GOTO 220 IF (J.NE.JNUL) GOTO 210 IF (JREST.EQ.1) WRITE(LUB,8090)(JDRAW(K),K=1,JXSCA) IF (JREST.NE.1) WRITE(LUB,8100)(JDRAW(K),K=1,JXSCA) GOTO 240 210 IF (J.NE.J2M.OR.(JPL2.EQ.0.AND.YKLE.NE.YMAX)) GOTO 220 IF (JREST.EQ.1) WRITE(LUB,8110) (JDRAW(K),K=1,JXSCA) IF (JREST.NE.1) WRITE(LUB,8120) (JDRAW(K),K=1,JXSCA) GOTO 240 220 IF(J.NE.JYSCA) GOTO 230 WRITE(LUB,8130) YKLE2,(JDRAW(K),K=1,JXSCA) GOTO 240 230 IF (JREST.EQ.1) WRITE(LUB,8140)(JDRAW(K),K=1,JXSCA) IF (JREST.NE.1) WRITE(LUB,8150)(JDRAW(K),K=1,JXSCA) 240 CONTINUE WRITE(LUB,8075) IF (JPLT.NE.2.OR.JPLXY.EQ.1) WRITE(LUB,8160) XMIN,XMAX IF ((JPLT.EQ.1.OR.JPLT.EQ.3).AND.JPLXY.EQ.0) 1 WRITE(LUB,8170) LAB(NVAD) IF (JPLXY.EQ.1) WRITE(LUB,8180) LAB(1) IF (JPLT.NE.2.OR.JPLXY.EQ.1) GOTO 250 J1=XMIN JN=XMAX WRITE(LUB,8190) J1,JN WRITE(LUB,8200) 8000 FORMAT(//16X,A60) 8010 FORMAT(//16X,A60//29X,32H--- L E A S T S Q U A R E S ---) 8020 FORMAT(//16X,A60// 1 20X,'--- L E A S T M E D I A N O F S Q U A R E S ---') 8025 FORMAT(//16X,A60// 1 22X,'--- L E A S T T R I M M E D S Q U A R E S ---') 8030 FORMAT(//16X,A60//29X,'--- R E W E I G H T E D',3X,7HL S ---) 8050 FORMAT(/' residual/scale',1X,3HI-+,10(5H----+),2H-I) 8060 FORMAT(/7X,8hresidual,3X,3HI-+,10(5H----+),2H-I) 8070 FORMAT(/6X,' observed'/6X,A10,2X,3HI-+,10(5H----+),2H-I) 8075 FORMAT(18X,1HI,53X,1HI) 8080 FORMAT(14X,3H2.5,1X,2H++,51A1,2H++) 8085 FORMAT(14X,3H2.5,1X,2HI+,51A1,2H+I) 8090 FORMAT(14X,3H0.0,1X,2H+-,51A1,2H-+) 8100 FORMAT(14X,3H0.0,1X,2HI-,51A1,2H-I) 8110 FORMAT(13X,4H-2.5,1X,2H++,51A1,2H++) 8120 FORMAT(13X,4H-2.5,1X,2HI+,51A1,2H+I) 8130 FORMAT(6X,E11.4,1X,1H+,1X,51A1,1X,1H+) 8140 FORMAT(18X,1H+,1X,51A1,1X,1H+) 8150 FORMAT(18X,1HI,1X,51A1,1X,1HI) 8160 FORMAT(18X,3HI-+,10(5H----+),2H-I/10X,E11.4,46X,E11.4) 8170 FORMAT(/49X,11h estimated ,A10//) 8180 FORMAT(/54X,9hobserved ,A10//) 8190 FORMAT(18X,3HI-+,10(5H----+),2H-I/17X,I4,46X,I4) 8200 FORMAT(/48X,25h index of the observation//) 250 RETURN END CC ----------------------------------------------------------------- CC *JQUANT* : ASKS FOR THE QUANTILE TO BE MINIMIZED CC ----------------------------------------------------------------- SUBROUTINE JQUANT(NCAS,NVAR,MROB,JQU,FACTOR,FACLMS,FACLTS) DIMENSION FACLMS(11),FACLTS(11) CHARACTER YN CC-----OPTIMAL BREAKDOWN VALUE IF H=(NCAS+NVAR+1)/2 JDEFAUL=(NCAS+NVAR+1)/2 JQU=JDEFAUL IF (MROB.EQ.1) THEN WRITE(*,8000) JDEFAUL ELSE WRITE(*,8010) JDEFAUL ENDIF JMIN=(NCAS/2)+1 JMAX=MAX((3*NCAS/4)+(NVAR+1)/4,JDEFAUL) IF (JMIN.NE.JMAX) THEN WRITE(*,8020) JDEFAUL 10 READ(*,8100) YN IF (YN.EQ.'y') YN='Y' IF (YN.EQ.'n') YN='N' IF (YN.NE.'Y'.AND.YN.NE.'N') THEN WRITE(*,8030) GOTO 10 ENDIF ENDIF IF (YN.EQ.'N') THEN JBRMIN=NBREAK(JMIN,NCAS,NVAR) JBRMAX=NBREAK(JMAX,NCAS,NVAR) WRITE(*,8040) JMIN,JBRMIN,JMAX,JBRMAX 20 READ(*,8110) JQU IF (JQU.LT.JMIN.OR.JQU.GT.JMAX) THEN WRITE(*,8030) GOTO 20 ENDIF ENDIF NQUANT=MIN(NINT(REAL(((JQU*1.0/NCAS)-0.5)*40))+1,11) IF (MROB.EQ.1) FACTOR=FACLMS(NQUANT) IF (MROB.EQ.2) FACTOR=FACLTS(NQUANT) 8000 FORMAT(//' The program will minimize the', I4,'th ordered', 1 ' squared residual.') 8010 FORMAT(//' The program will minimize the sum of the',I4, 1 ' smallest squared residuals.') 8020 FORMAT(' Do you agree with the default choice of ',I4,' ?', 1 ' Yes or no : ' ,$) 8030 FORMAT(/' Not allowed ! Enter your choice again : ',$) 8040 FORMAT(' Give a number between ',I4,' (then the breakdown', 1 ' value is ',I2,'%)'/ 1 ' and ',I4, 1 ' (then the breakdown value is ',I2,'%) : ',$) 8100 FORMAT(A1) 8110 FORMAT(I4) RETURN END CC ------------------------------------------------------------------- CC *NBREAK* : COMPUTES BREAKDOWN VALUE OF LMS/LTS CC ------------------------------------------------------------------- FUNCTION NBREAK(NQU,NCAS,NVAR) IF (NQU.LT.(NCAS+NVAR+1)/2) THEN NBREAK=(NQU-NVAR+1)*100/NCAS ELSE NBREAK=(NCAS-NQU+1)*100/NCAS ENDIF RETURN END