//ONTARIO JOB /*JOBPARM S=15,R=4096,L=16 /*SERVICE -4 //STEP1 EXEC FORTVCLG,RC=2048K,RL=2048K,RG=2048K,LIB=LIBD, // PARM.FORT='LANGLVL(77),NOMAP,NOXREF,OPTIMIZE(1)', // PARM.LKED='LET' //FORT.SYSIN DD * IMPLICIT REAL*8(A-H,O-Z) CHARACTER LABEL*35(6), OPTION*2(3) C REAL * 4 XBOX(2000) INTEGER KBOX, NIBOX(50), MAXL, IERBOX C DIMENSION TABLE(2000,4) DIMENSION DX(2000), DY(2000), DZ(2000), DT(2000) DIMENSION DL1(2000), DL2(2000), DL3(2000), DL4(2000) C DIMENSION TZERO(3) C DIMENSION A(2000,3), T(3), WORK(2000) DIMENSION S(3,3), ST(3,3), B(3,3), C(3,3) DIMENSION SINV(3,3), IW1(3), IW2(3) DIMENSION V1(2000), ARRAY(2000), ARRAY1(2000) DIMENSION COV(3,3) , XI(3,1), DXI(3,1) C DIMENSION V(3,3), Y(2000,3), YNORM(2000) DIMENSION W(3,3) DIMENSION ALFA(3,1), BETA(1,3), OUT(3,3) DIMENSION COVUP(3,3), H(3), VH(3) C DIMENSION COVCAM(3,3), INDEX(2000) DIMENSION VHUBE(3,3), COVHUB(3,3) DIMENSION VCAM(3,3), COVMVT(3,3), TCAM(3) DIMENSION VTRM(3,3), AA(2000,3), INDX(2000) DIMENSION VMAR(3,3), COVMAR(3,3) C DIMENSION XELL(100), YELL(100) DIMENSION XEL1(100), YEL1(100) DIMENSION D(3), EGV(3,3), WK(18) DIMENSION D2(2), EGV2(2,2) DIMENSION VECS(6) C DIMENSION VHUB(200), VCBL1(200), VCBL2(200) DIMENSION VCLS1(200), VCLS2(200), VCLS3(200) DIMENSION VHUB1(200), VHUB2(200), VHUB3(200) DIMENSION VCBL1A(200), VCBL1B(200), VCBL1C(200) DIMENSION VCBL2A(200), VCBL2B(200), VCBL2C(200) DIMENSION VMVT1(200), VMVT2(200), VMVT3(200) DIMENSION VMLC1(200), VMLC2(200), VMLC3(200) DIMENSION VMVT(200), VMLC(200), VCLS(200) DIMENSION SCALE(6,200), VARLOG(6), EFF(6) DIMENSION DIREC(3), DIR(3), DIRCOS(200) C LABEL(1)= 'CLASSICAL' LABEL(2)= ' HUBER ' LABEL(3)= 'CAMPBELL1' LABEL(4)= 'CAMPBELL2' LABEL(5)= 'M V T' LABEL(6)= 'MARONNA ' C OPTION(1) = ' A' OPTION(2) = ' B' OPTION(3) = ' C' C IP = 0 IPRINT = 6 ITEMP = 15 IDISK = 22 IREAD = 5 NPARA = 3 IRDA = 2000 PI = 3.14159D0 MAXL = 129 C---------------------------------------------------------- C IOPT = 1 RAW DATA C IOPT = 2 OUTLIERS HAVE BEEN REMOVED C IOPT = 3 GOOD GEOMETRICAL CONDITIONS N=400 ONLY C IEFF = 0 IF NO EFFICIENCY COMPUTATIONS TAKE PLACE C IEFF > 0 IF EFFICIENCY COMPUTATIONS TAKE PLACE C--------------------------------------------------------- IOPT = 1 IEFF = 0 TZERO(1) = 0.D0 TZERO(2) = 0.D0 TZERO(3) = 0.D0 IGEN = 1 ISCAL = 6 ISAMPL= 20.D0 RANGE = 1167.0D0 C------------------------------------------------------------------ C READ DATA C----------------------------------------------------------------- OPEN(IDISK,ERR=3000,STATUS='OLD',ACCESS='SEQUENTIAL', # IOSTAT=IERROR) C 9 DO 1 K=1,2000 READ(IDISK,142,ERR=2010,IOSTAT=IERROR,END=2) INDEX(K), DX(K), # DY(K), DZ(K), DT(K), DL1(K), DL2(K), DL3(K), DL4(K) 1 CONTINUE 2 NDATA = K - 1 NOBS = NDATA C--------------------------------------------------- C FORM DATA TABLE C SUBTRACT MEDIANS C-------------------------------------------------- DO 3 I= 1,NDATA TABLE(I,1) = DX(I) - 1.8770D0 TABLE(I,2) = DY(I) - 2.2800D0 TABLE(I,3) = DZ(I) - 0.636D0 RT=DSQRT( TABLE(I,1)**2 +TABLE(I,2)**2+TABLE(I,3)**2) TABLE(I,4) = ( RT) 3 CONTINUE IF(IOPT .EQ. 2) THEN C-----------------------------WATCH OUT ------------ CRIT = 0.06D0 C------------------------------------- MODE = 2 B1 = 2.0D0 ACBL = R0CAM(NPARA, B1) BCBL = 1.25D0 EPSI= 10.D-3 CALL COVMTX(TABLE, IRDA, NDATA, NPARA, T, COVCAM) CALL COVML(TABLE, IRDA, NDATA, NPARA, T, ACBL, BCBL, EPSI, # VCAM, COVCAM ,IP, MODE) C-------------------------------------------------------------- C REJECT OUTLIERS C-------------------------------------------------------------- CALL GROSS(TABLE, IRDA, NDATA, NPARA, ACBL, BCBL, CRIT, T, # VCAM, COVCAM, INDEX, ICLEAN) NDATA = ICLEAN NOBS = NDATA RANGE = DFLOAT(ICLEAN) ENDIF C IF( IOPT .EQ. 3) THEN RANGE = 400.D0 NDATA = 400 NOBS = NDATA ENDIF WRITE(6,170) OPTION(IOPT), NDATA, ISAMPL, IGEN C******************************************************************** DO 9999 K=1, IGEN CALL GENRDM(TABLE, IRDA, NPARA, 4, ISAMPL, RANGE) DO 11 I=1, ISAMPL C READ(5,*) TABLE(I,1), TABLE(I,2), TABLE(I,3) 11 WRITE(6,111) I, TABLE(I,1), TABLE(I,2), TABLE(I,3) 111 FORMAT(2X,I3,3(F10.2)) NDATA = ISAMPL C------------------------------------------------------------------- C START WITH CLASSICAL ESTIMATES FOR LOCATION C-------------------------------------------------------------------- 77 ITER = 0 WRITE(6,150) LABEL(1), K CALL LOCATE(TABLE, IRDA, NDATA, NPARA, T) CALL COVMTX(TABLE, IRDA, NDATA, NPARA, T, COVHUB) CALL MOUTD(COVHUB, 3, 3, 3) CALL STORE(COVHUB, 3, VECS) CALL EIGRS(VECS, 3, 1, D, EGV, 3, WK, IERROR) CALL MOUTD(D, 3, 3, 1) CALL MOUTD(EGV, 3, 3, 3) C CALL MOUTD(T, 3, 3, 1) VCLS1(K) = DLOG( D(1) ) VCLS2(K) = DLOG( D(2) ) VCLS3(K) = DLOG( D(3) ) C VCLS(K) = VCLS1(K) + VCLS2(K) + VCLS3(K) C C------------- ASYMPTOTIC EFFICIENCY ------------------- IF( IEFF. GT. 0) THEN CALL CHOLSK(COVHUB, 3, 3, S, NPARA) C CALL MOUTD(S, 3, 3, 3) CALL XNVLOW(S, NPARA, NPARA, V) C CALL MOUTD(V, 3, 3, 3) CALL YOUT(TABLE, IRDA, NDATA, NPARA, V, TZERO, Y, YNORM) DO 33 I=1, NDATA R = YNORM(I) ARRAY(I) = R** 2 ARRAY1(I) = 2.D0 * R 33 CONTINUE CALL ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) WRITE(6,888) ALAMDA 888 FORMAT(3X,'LAMDA =',F12.6,2X,'CLASSICAL') ENDIF C------------------2D CASE---------------------------- VECS(1) = COVHUB(1,1) VECS(2) = COVHUB(2,1) VECS(3) = COVHUB(2,2) C CALL EIGRS(VECS, 3, 1, D, EGV, 2, WK, IERROR) C CALL MOUTD(D,3,3,1) C CALL MOUTD(EGV,3,3,3) IF(K .EQ. 1) THEN DIREC(1) = EGV(1,2) DIREC(2) = EGV(2,2) DIREC(3) = EGV(3,2) ENDIF DIR(1) = EGV(1,2) DIR(2) = EGV(2,2) DIR(3) = EGV(3,2) CALL INNER(DIREC, DIR, NPARA, DIRCOS(K)) DIRCOS(K) = DABS(DIRCOS(K)) C C CALL MOUTD(D, 3, 3, 1) C CALL MOUTD(EGV, 3, 3, 3) C CALL MOUTD(COVHUB, 3, 3, 3) C C------------------------------------------------------------------- C HUBER'S APPROACH C------------------------------------------------------------------- WRITE(6,150) LABEL(2), K ITER = 0 5 MODE = 1 AH = 0.D0 BH = 2.50D0 CALL TAF2(AH,BH,2) EPSI= 10.D-3 CALL COVML(TABLE, IRDA, NDATA, NPARA, T, AH, BH, EPSI, # VHUBE, COVHUB,IP, MODE) C---------------------------------LOCATION ---------------------- AM = 0.D0 BM = 2.0D0 CALL LOCML(TABLE, IRDA, NDATA, NPARA,VHUBE, AM, BM, # EPSI, T, H, VHNRM, MODE,IP) 91 ITER = ITER + 1 IF(ITER. EQ. 1) GO TO 5 CALL MOUTD(T,3,3,1) C-------------------- ASYMPTOTIC EFFICIENCY ------------ IF( IEFF .GT. 0) THEN CALL YOUT(TABLE, IRDA, NDATA, NPARA, VHUBE, TZERO, Y, YNORM) DO 331 I=1,NDATA R = YNORM(I) ARRAY(I) = UBR( AH, BH, R) ARRAY1(I)= UBRDRV(AH, BH, R) 331 CONTINUE CALL ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) WRITE(6,881) ALAMDA 881 FORMAT(4X,'LAMDA =', F10.6,' HUBER ') ENDIF CALL STORE(COVHUB, 3, VECS) CALL EIGRS(VECS, 3, 1, D, EGV, 3, WK, IERROR) VHUB1(K) = DLOG( D(1) ) VHUB2(K) = DLOG( D(2) ) VHUB3(K) = DLOG( D(3) ) VHUB(K) = VHUB1(K) + VHUB2(K) + VHUB3(K) C VECS(1) = COVHUB(1,1) VECS(2) = COVHUB(3,1) VECS(3) = COVHUB(3,3) CALL MOUTD(D, 3, 3, 1) CALL MOUTD(EGV, 3, 3, 3) CALL MOUTD(COVHUB, 3, 3, 3) IF(1 .EQ. 1) GO TO 2222 C------------------------------------------------------------------ C CAMPBELL'S APPROACH C------------------------------------------------------------------ C WRITE(6,150) LABEL(3), K MODE = 2 B1 = 2.0D0 ACBL = R0CAM(NPARA, B1) BCBL = 1.25D0 EPSI= 10.D-3 ITER = 0 CALL SETD(COVCAM, NPARA,COVHUB) 6 CALL COVML(TABLE, IRDA, NDATA, NPARA, T, ACBL, BCBL, EPSI, # VCAM, COVCAM ,IP, MODE) C------------------------------------------------------------- C LOCATION C------------------------------------------------------------- CALL LOCML(TABLE, IRDA, NDATA, NPARA,VCAM, ACBL, BCBL, # EPSI, T, H, VHNRM, MODE,IP) ITER = ITER + 1 IF(ITER .EQ. 1) GO TO 6 CALL STORE(COVCAM, 3, VECS) CALL EIGRS(VECS, 3, 1, D, EGV, 3, WK, IERROR) VCBL1(K) = DLOG( D(1) ) + DLOG( D(2) ) + DLOG( D(3) ) C--------------------- ASYMPTOTIC EFFICIENCY ---------- IF( IEFF .GT. 0) THEN CALL YOUT(TABLE, IRDA, NDATA, NPARA, VCAM, TZERO, Y, YNORM) DO 332 I=1, NDATA R = YNORM(I) UCBL = WSCAM(ACBL, BCBL, R) ARRAY(I) = UCBL ** 2 ARRAY1(I) = CAMDRV(ACBL, BCBL, R) 332 CONTINUE CALL ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) WRITE(6,882) ALAMDA 882 FORMAT(4X,'LAMDA =', F10.6,' CAMPBELL FIRST ') ENDIF C CALL MOUTD(D, 3, 3, 1) C CALL MOUTD(EGV, 3, 3, 3) VCBL1A(K) = DLOG( D(1) ) VCBL1B(K) = DLOG( D(2) ) VCBL1C(K) = DLOG( D(3) ) C C------------------------------------------------------------ C CAMPBELL'S SECOND APPROACH C------------------------------------------------------------ C WRITE(6,150) LABEL(4), K BCBL = 50000.D0 ITER = 0 CALL SETD(COVCAM, NPARA,COVHUB) 61 CALL COVML(TABLE, IRDA, NDATA, NPARA, T, ACBL, BCBL, EPSI, # VCAM, COVCAM ,IP, MODE) CALL LOCML(TABLE, IRDA, NDATA, NPARA,VCAM, ACBL, BCBL, # EPSI, T, H, VHNRM, MODE,IP) ITER = ITER + 1 IF(ITER .EQ. 1) GO TO 61 CALL STORE(COVCAM, 3, VECS) CALL EIGRS(VECS, 3, 1, D, EGV, 3, WK, IERROR) VCBL2(K) = DLOG( D(1) ) + DLOG( D(2) ) + DLOG( D(3) ) C C--------------------- ASYMPTOTIC EFFICIENCY ---------- IF( IEFF .GT. 0) THEN CALL YOUT(TABLE, IRDA, NDATA, NPARA, VCAM, TZERO, Y, YNORM) DO 333 I=1, NDATA R = YNORM(I) UCBL = WSCAM(ACBL, BCBL, R) ARRAY(I) = UCBL ** 2 ARRAY1(I) = CAMDRV(ACBL, BCBL, R) 333 CONTINUE CALL ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) WRITE(6,883) ALAMDA 883 FORMAT(4X,'LAMDA =', F10.6,' CAMPBELL SECOND ') ENDIF C C CALL MOUTD(D, 3, 3, 1) C CALL MOUTD(EGV, 3, 3, 3) VCBL2A(K) = DLOG( D(1) ) VCBL2B(K) = DLOG( D(2) ) VCBL2C(K) = DLOG( D(3) ) C C------------------------------------------------------------ C MULTIVARIATE TRIMMING C------------------------------------------------------------ C WRITE(6,150) LABEL(5), K ALF = 0.10D0 CALL TRIM(TABLE, IRDA,NDATA, NPARA, T, ALF, COVMVT, AA,INDX) C C--------------------- ASYMPTOTIC EFFICIENCY ------------- IF( IEFF .GT. 0) THEN CALL CHOLSK(COVMVT, 3, 3, S, NPARA) CALL XNVLOW(S, 3, 3, VTRM ) C KTRIM = ALF * NDATA + 0.5D0 ITRIM = NDATA - KTRIM CALL YOUT(AA , IRDA, ITRIM, NPARA, VTRM, TZERO, Y, YNORM) DO 50 I=1,NDATA ARRAY(I) = 0.D0 ARRAY1(I)= 0.D0 IF(I .GT. ITRIM) YNORM(I) =0.D0 50 CONTINUE DO 334 I=1, ITRIM R = YNORM(I) ARRAY(I) = R ** 2 ARRAY1(I) = 2 * R 334 CONTINUE CALL ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) WRITE(6,884) ALAMDA, ITRIM 884 FORMAT(4X,'LAMDA =', F10.6,' M V T ',3X,'ITRIM=',I8) ENDIF C CALL STORE(COVMVT, 3, VECS) CALL EIGRS(VECS, 3, 1, D, EGV, 3, WK, IERROR) C CALL MOUTD(D, 3, 3, 1) C CALL MOUTD(EGV, 3, 3, 3) VMVT1(K) = DLOG( D(1) ) VMVT2(K) = DLOG( D(2) ) VMVT3(K) = DLOG( D(3) ) C VMVT(K) = DLOG( D(1) ) + DLOG( D(2) ) + DLOG( D(3) ) C------------------------------------------------------------ C MARONNA'S APPROACH C------------------------------------------------------------ C WRITE(6,150) LABEL(6), K MODE = 3 AMAR = 0.D0 BMAR = 0.D0 ITER = 0 CALL COVMTX(TABLE, IRDA, NDATA, NPARA, T, COVMAR) 7 CALL COVML(TABLE, IRDA, NDATA, NPARA, T, AMAR, BMAR, EPSI, # VMAR, COVMAR ,IP, MODE) C----------------------------LOCATION------------------------- CALL LOCML(TABLE, IRDA, NDATA, NPARA,VMAR, AMAR, BMAR, # EPSI, T, H, VHNRM, MODE,IP) ITER = ITER + 1 IF(ITER .EQ. 1) GO TO 7 C C--------------------- ASYMPTOTIC EFFICIENCY ---------- IF(IEFF .GT. 0) THEN CALL YOUT(TABLE, IRDA, NDATA, NPARA, VMAR, TZERO, Y, YNORM) DO 335 I=1, NDATA R = YNORM(I) UMAR = 4.D0 / ( 1.D0 + (1.D0 / (R**2) ) ) ARRAY(I) = UMAR ARIT = 8.D0 * R PARO = (1.D0 + R**2) ** 2 ARRAY1(I) = ARIT / PARO 335 CONTINUE CALL ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) WRITE(6,885) ALAMDA 885 FORMAT(4X,'LAMDA =', F10.6,' MARONNA ') ENDIF CALL STORE(COVMAR, 3, VECS) CALL EIGRS(VECS, 3, 1, D, EGV, 3, WK, IERROR) C CALL MOUTD(D, 3, 3, 1) C CALL MOUTD(EGV, 3, 3, 3) VMLC1(K) = DLOG( D(1) ) VMLC2(K) = DLOG( D(2) ) VMLC3(K) = DLOG( D(3) ) VMLC(K) = DLOG( D(1) ) + DLOG( D(2) ) + DLOG( D(3) ) C NDATA = NOBS 9999 CONTINUE C*************************************************************** WRITE(6,170) OPTION(IOPT), NDATA, ISAMPL, IGEN IF(IEFF. GT. 0) GO TO 2020 WRITE(6,180) C------------------------------------- CALL SORT8(VCLS, IGEN ) CALL SORT8(VHUB, IGEN ) CALL SORT8(VCBL1,IGEN ) CALL SORT8(VCBL2,IGEN ) CALL SORT8(VMVT, IGEN ) CALL SORT8(VMLC, IGEN ) DO 999 I=1, IGEN WRITE(15,*) VCLS(I),VHUB(I), VCBL1(I), VCBL2(I),VMVT(I), # VMLC(I) WRITE(6,160) I, VCLS(I),VHUB(I), VCBL1(I), VCBL2(I),VMVT(I), # VMLC(I), DIRCOS(I) 999 CONTINUE C-------------------------------------------------------------- C COMPUTE EFFICIENCIES C-------------------------------------------------------------- CALL GENMTX(VCLS, VHUB,VCBL1, VCBL2, VMVT, VMLC, # SCALE , ISCAL, IGEN) CALL VRNLOG(SCALE, ISCAL, IGEN, VARLOG, EFF) WRITE(6,200) DO 80 I=1, ISCAL WRITE(6,190) I, LABEL(I), VARLOG(I), EFF(I) 80 CONTINUE CALL STRDEV(VCLS1, IGEN, SCL1) CALL STRDEV(VCLS2, IGEN, SCL2) CALL STRDEV(VCLS3, IGEN, SCL3) C CALL STRDEV(VHUB1, IGEN, SHU1) CALL STRDEV(VHUB2, IGEN, SHU2) CALL STRDEV(VHUB3, IGEN, SHU3) C CALL STRDEV(VCBL1A, IGEN, SCB1) CALL STRDEV(VCBL1B, IGEN, SCB2) CALL STRDEV(VCBL1C, IGEN, SCB3) C CALL STRDEV(VCBL2A, IGEN, SC1B) CALL STRDEV(VCBL2B, IGEN, SC2B) CALL STRDEV(VCBL2C, IGEN, SC3B) C CALL STRDEV(VMVT1, IGEN, SMV1) CALL STRDEV(VMVT2, IGEN, SMV2) CALL STRDEV(VMVT3, IGEN, SMV3) C CALL STRDEV(VMLC1, IGEN, SML1) CALL STRDEV(VMLC2, IGEN, SML2) CALL STRDEV(VMLC3, IGEN, SML3) C WRITE(6,220) LABEL(1), SCL1**2, SCL2**2, SCL3**2 WRITE(6,220) LABEL(2), SHU1**2, SHU2**2, SHU3**2 WRITE(6,220) LABEL(3), SCB1**2, SCB2**2, SCB3**2 WRITE(6,220) LABEL(4), SC1B**2, SC2B**2, SC3B**2 WRITE(6,220) LABEL(5), SMV1**2, SMV2**2, SMV3**2 WRITE(6,220) LABEL(6), SML1**2, SML2**2, SML3**2 C------------------------------------------------------------------ C CLOSE FILE ON DISK C------------------------------------------------- GO TO 2020 3000 WRITE(IPRINT,9040) IERROR GO TO 2020 2010 WRITE(IPRINT,9030) IERROR GO TO 2020 2020 CLOSE(IDISK,ERR=4010,STATUS='KEEP',IOSTAT=IERROR) GO TO 2222 C--------------------------------------------------------------------- 63 FORMAT(2X,'ALFA=', F10.3,' BETA=',F10.3,' TAF=', F10.5, # ' ETA=',F10.3,' ZETA=',F10.3,' ANGLE=',F10.5) 142 FORMAT(I8,1X,3(F10.3,2X),F20.4,F20.4,1X,3(F14.2)) 150 FORMAT('1',3X,A35,3X,'SAMPLE NUMBER=',I6,//) 160 FORMAT(' ',3X,I5, 6(F10.3,2X), F12.6) 170 FORMAT('1','LOGARITHM OF ELLIPSOIDAL VOLUMES IN SITUATION=',A2, # 2X,'NDATA=',I6,3X,'SAMPLE SIZE=',I6,3X,'IGEN=',I6) 180 FORMAT(2X,'NUMBER',3X,'CLASSICAL',5X,'HUBER',6X,'CAMPBELL1',4X # ,'CAMPBELL2',4X, 'TRIMMING',5X,'MARONNA') 190 FORMAT(2X,I5,1X,A35,F10.5,3X,F7.3,//) 200 FORMAT('1','NUMBER',' COVARIANCE ESTIMATE VAR OF LOG(VOL #UME) EFFICIENCY',//) 220 FORMAT(2X,A20,3(F10.4,2X)) 300 FORMAT('1','+++++++++++++++++++++++++++++++++++++++++++++++++++++ #+++ CAMPBELL APPROACH ++++++++++') 9030 FORMAT(' ','ERROR',I4,3X,'WHILE READING FILE OF RESULTS') 9040 FORMAT(' ','ERROR',I4,3X,'WHILE OPENING FILE OF RESULTS') 9070 FORMAT(' ','ERROR',I4,3X,'WHILE CLOSING FILE OF RESULTS') C------------------------------------------------------ 4010 WRITE(IPRINT,9070) IERROR 2222 CONTINUE STOP END C FUNCTION UBRDRV(AH, BH, R) IMPLICIT REAL * 8(A-H, O-Z) C------------------------------------------------ C DERIVATIVE OF HUBER'S WEIGHT FOR COVARIANCE C AUTHOR : S. MERTIKAS C------------------------------------------------- IF( R .GE. AH .AND. R .LE. BH) THEN UBRDRV = 2.D0 * R RETURN ENDIF UBRDRV = 0.D0 RETURN END C SUBROUTINE ASYMPT(ARRAY, ARRAY1, NDATA, NPARA, YNORM, ALAMDA) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION ARRAY(NDATA), ARRAY1(NDATA), YNORM(NDATA) C------------------------------------------------------------------- C COMPUTES THE LAMDA COEEFICIENT IN THE ASYMPTOTIC EFFICIENCY C OF COVARIANCE MATRICES (HUBER,1977) C C ARRAY = VECTOR CONTAINING WEIGHTS FOR M-ESTIMATION IN COVARIANCE C ARRAY1= VECTOR CONTAINING DERIVATIVES OF PREVIOUS WEIGHTS C NPARA = NUMBER OF PARAMETERS C NDATA = NUMBER OF DATA C YNORM = VECTOR CONTAINING THE SQUARE-ROOTS OF DISTANCE C C ALAMDA = LAMDA COEFFICIENT (OUTPUT) C C AUTHOR : S. MERTIKAS C DATE : MARCH,1987 C-------------------------------------------------------------------- SUMA = 0.D0 SUMB = 0.D0 EN = DFLOAT(NDATA) C---------------COMPUTE NUMERATOR------------------------ DO 10 I=1, NDATA R2 = ARRAY(I) ** 2 SUMA = SUMA + R2 10 CONTINUE C--------------COMPUTE DENOMINATOR----------------------- DO 20 I=1,NDATA Q = ARRAY1(I) / (DFLOAT(NPARA)) R = YNORM(I) SUMB = SUMB + Q * R + ARRAY(I) 20 CONTINUE DEN = SUMB ** 2 C WRITE(6,6) SUMA, DEN, EN 6 FORMAT(5X,'SUM NUMERATOR=',F12.3,' SUM DENOM=',F15.2, # 'NDATA=',F8.1) ALAMDA = (SUMA / DEN) * EN RETURN END C SUBROUTINE MINVD (A,IA,MA,DETA,IR,IC) REAL*8A(IA,MA),DETA,PIV,PIV1,TEMP,TEST,X,DABS DIMENSIONIR(MA),IC(MA) C A CONTAINS MATRIX TO BE INVERTED C IA IS THE ROW DIENSION OF THE MATRIX A C MA IS THE SIZE OF THE MATRIX TO BE INVERTED. C DETA IS THE DETERMINANT C IR,IC ARE 2 WORK VECTORS. DO1I=1,MA IR(I)=0 1 IC(I)=0 DETA=1.0D0 DO123IJKL=1,MA I=0 J=0 TEST=0.0D0 DO55K=1,MA IF(IR(K).NE.0)GOTO55 DO4L=1,MA IF(IC(L).NE.0)GOTO4 X=DABS(A(K,L)) IF(X.LT.TEST)GOTO4 I=K J=L TEST=X 4 CONTINUE 55 CONTINUE PIV=A(I,J) DETA=PIV*DETA IF(PIV.EQ.0.D0)GOTO17 IR(I)=J IC(J)=I PIV=1.D0/PIV DO5K=1,MA 5 A(I,K)=A(I,K)*PIV A(I,J)=PIV DO9K=1,MA IF(K.EQ.I)GOTO9 PIV1=A(K,J) 6 DO8L=1,MA 8 A(K,L)=A(K,L)-PIV1*A(I,L) A(K,J)=PIV1 9 CONTINUE PIV1=A(I,J) DO11K=1,MA 11 A(K,J)=-PIV*A(K,J) A(I,J)=PIV1 123 CONTINUE 12 DO16I=1,MA K=IC(I) M=IR(I) IF(K.EQ.I)GOTO16 DETA=-DETA DO14L=1,MA TEMP=A(K,L) A(K,L)=A(I,L) 14 A(I,L)=TEMP DO15L=1,MA TEMP=A(L,M) A(L,M)=A(L,I) 15 A(L,I)=TEMP IC(M)=K IR(K)=M 16 CONTINUE RETURN 17 DETA=0.0D0 RETURN END C SUBROUTINE ROTATE(S1, S2, S12, ANGLE) IMPLICIT REAL * 8(A-H,O-Z) PI = 3.14159D0 DY = 2.D0 * S12 DX = S1 - S2 ANGLE = 0.5D0* DATAN2(DY, DX) RETURN END C SUBROUTINE INNER(ARRAY, ARRAY1, NPARA, COSA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ARRAY(NPARA), ARRAY1(NPARA) C------------------------------------------------ C COMPUTES THE INNER PRODUCT OF TWO VECTORS C-------------------------------------------------- COSA = 0.D0 DO 10 I=1, NPARA COSA = COSA + ARRAY(I) * ARRAY1(I) 10 CONTINUE RETURN END C SUBROUTINE XMEAN(ARRAY,NDATA,XBAR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY(NDATA) XBAR = 0.D0 DO 1 I=1,NDATA 1 XBAR=XBAR+ARRAY(I) XBAR=XBAR/DFLOAT(NDATA) RETURN END C SUBROUTINE STRDEV(ARRAY, NDATA, SD) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION ARRAY( NDATA ) C----------------------------------------------------------------- C COMPUTES STANDARD DEVIATION OF A DATA SET ARRAY C C INPUT : ARRAY = VECTOR OF SAMPLE DATA C NDATA = SAMPLE SIZE C OUTPUT: SD = STANDARD DEVIATION C C AUTHOR: S. MERTIKAS C DATE : AUGUST,1986 C--------------------------------------------------------------- CALL XMEAN(ARRAY,NDATA,XBAR) SUM = 0.0D0 DO 10 I =1, NDATA SUM = SUM + (ARRAY(I) - XBAR) ** 2 10 CONTINUE SD = SUM / DFLOAT(NDATA - 1) SD = DSQRT( SD ) RETURN END SUBROUTINE SETV(ARRAY, N, IFIRST, ILAST, ARRAY1, NDATA) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY( N ), ARRAY1(NDATA) C J = 0 DO 1 I = IFIRST, ILAST J = J + 1 ARRAY( I ) = ARRAY1( J ) 1 CONTINUE RETURN END C SUBROUTINE VRNLOG(SCALE, ISCAL, IGEN, VARLOG, EFF) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION SCALE(ISCAL, IGEN), VARLOG(ISCAL), SLOG(200) DIMENSION YY(200), EFF(ISCAL) C------------------------------------------------------------------- C COMPUTES THE VARIANCE OF THE LOGARITHM OF VOLUME C A CONFIDENCE ELLIPSOID C AND CORRESPONDING EFFICIENCIES C AUTHOR : S. MERTIKAS C DATE : AUGUST,1986 C------------------------------------------------------------------- DO 10 I=1, ISCAL DO 20 J=1, IGEN SLOG(J) = SCALE(I,J) 20 CONTINUE CALL STRDEV( SLOG, IGEN, SIGMA ) VARLOG(I) = SIGMA ** 2 10 CONTINUE CALL SETV(YY, ISCAL,1,ISCAL, VARLOG, ISCAL) CALL SORT8(YY, ISCAL) VMIN = YY(1) DO 30 I=1, ISCAL EFF(I) = (VMIN*100.D0) / VARLOG(I) 30 CONTINUE RETURN END C SUBROUTINE GENMTX(VCLS, VHUB,VCBL1, VCBL2, VMVT, VMLC, # SCALE , ISCAL, IGEN) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION VHUB(IGEN), VCBL1(IGEN), VCBL2(IGEN) DIMENSION VMVT(IGEN), VMLC(IGEN), VCLS(IGEN) DIMENSION SCALE(ISCAL, IGEN) C--------------------------------------------------------------------- C GENERATES MATRIX OF ALL SCALE ESTIMATES FOR COMPUTING EFFICIENCY C--------------------------------------------------------------------- DO 10 K=1, IGEN SCALE(1,K) = VCLS(K) SCALE(2,K) = VHUB(K) SCALE(3,K) = VCBL1(K) SCALE(4,K) = VCBL2(K) SCALE(5,K) = VMVT(K) SCALE(6,K) = VMLC(K) 10 CONTINUE RETURN END C SUBROUTINE SETD(COVCAM, NPARA, COVHUB) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION COVCAM(NPARA,NPARA), COVHUB(NPARA, NPARA) DO 10 I= 1, NPARA DO 20 J=1,NPARA COVCAM(I,J) = COVHUB(I,J) 20 CONTINUE 10 CONTINUE RETURN END C SUBROUTINE GENRDM(A, IRDA, NPARA, IGEN, ISAMPL, RANGE) IMPLICIT REAL * 8(A-H, O-Z) DIMENSION A( IRDA, NPARA ) C----------------------------------------------------------------- C GENERATES A MATRIX OF SIZE ISAMPL OF RANDOMLY C SELECTED VALUES FROM A(IRDA, NPARA) C---------------------------------------------------------------- IY = IGEN DO 10 I=1, ISAMPL RANDOM = URAND(IY) TEST = RANDOM * RANGE IR = TEST IF(IR .EQ. 0) IR = 1 A(I, 1) = A(IR, 1) A(I, 2) = A(IR, 2) A(I, 3) = A(IR, 3) C WRITE(6,*) I, IR, TEST, A(I,1) 10 CONTINUE RETURN END C SUBROUTINE STORE(COV, N, VECS) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION COV(N,N), VECS(6) K = 0 DO 10 I=1, N DO 20 J=1, I K = K + 1 VECS(K) = COV(I,J) 20 CONTINUE 10 CONTINUE RETURN END C SUBROUTINE TAF2(A, B, NPARA) IMPLICIT REAL * 8(A-H, O-Z) REAL SP, SP2 INTEGER IER T = 1.00D0 P = DFLOAT(NPARA) SP = SNGL(P) SP2= SP + 2.0 DO 10 I= 1, 400 T = T + 0.001D0 RA = (A / T) ** 2 RB = (B / T) ** 2 CALL MDCH(SNGL(RA), SP,DBLE(P1), IER) TRA = (A**2) * P1 CALL MDCH(SNGL(RB), SP,DBLE( P2), IER) TRB = (B**2) * (1.D0 - P2) CALL MDCH(SNGL(RB), SP2, DBLE(P3), IER) CALL MDCH(SNGL(RA), SP2, DBLE(P4), IER) TRC = (T ** 2) * P * ( P3- P4) RES = TRA + TRB + TRC WRITE(6,100) RES, T, I 100 FORMAT(4X,' COMP PARAMETER=',F10.3,' TAF=',F10.4,' I=',I4) 10 CONTINUE RETURN END C SUBROUTINE ELIPSE(XELL, YELL,NP, ETA, ZETA, AEL, BEL, ANGLE) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION XELL(NP), YELL(NP) C--------------------------------------------------------- C COMPUTE X AND Y COORDINATES FOR DRAWING AN ELLIPSE C C INPUT : N = NUMBER OF POINTS (<100) C : ETA = X-ORIGIN C : ZETA= Y-ORIGIN C : AEL = MAJOR-AXIS OF ELLIPSE C : BEL = MINOR-AXIS OF ELLIPSE C OUTPUT : XEL = VECTOR OF X-COORDINATES C : YEL = VECTOR OF Y-COORDINATES C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C--------------------------------------------------------- STEP = AEL / (22.D0) XI = 0.D0 DO 10 I =1, 12 DX = XI D = 1.D0 - (DX / AEL) ** 2 DY = BEL * DSQRT( D ) C XELL(I) = DX YELL(I) = DY C J = 51 -I XELL(J) = XELL(I) YELL(J) = - DY C L = I + 50 XELL(L) = - DX YELL(L) = - DY C K = 101 - I XELL(K) = - DX YELL(K) = YELL(I) C XI = XI + STEP 10 CONTINUE YI = (BEL / 2.D0) STEP = BEL / 24.D0 DO 20 I =13,25 DY = YI D = 1.D0 - (DY / BEL) ** 2 DX = AEL * DSQRT( D ) C XELL(I) = DX YELL(I) = DY C J = 51 -I XELL(J) = XELL(I) YELL(J) = - DY C L = I + 50 XELL(L) = - DX YELL(L) = - DY C K = 101 - I XELL(K) = - DX YELL(K) = YELL(I) C YI = YI - STEP 20 CONTINUE C------------------------------------ C ROTATE C------------------------------------ DO 30 I =1 , 100 XE = DCOS(ANGLE) * XELL(I) -DSIN(ANGLE) * YELL(I) YE = DSIN(ANGLE) * XELL(I) +DCOS(ANGLE) * YELL(I) WRITE(6,66) XE, YE, XELL(I), YELL(I), I 66 FORMAT(2X,'XNEW=',F10.3,' YNEW=',F10.3,' X=',F10.3,' Y=',F10.3, # ' I=',I5) XELL(I) = XE + ETA YELL(I) = YE + ZETA 30 CONTINUE RETURN END C FUNCTION WMAR(NPARA, F, R) IMPLICIT REAL * 8 (A-H,O-Z) C------------------------------------------------------ C COMPUTES MARONNA'S WEIGHT C------------------------------------------------------ W = (DFLOAT(NPARA) + F ) / ( F + R ** 2) WS = (DFLOAT(NPARA) + F ) / ( F ) WMAR= W / WS RETURN END C SUBROUTINE TRIM(A, IRDA,NDATA, NPARA, T, ALF, COV,AA,INDX) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION A(IRDA,NPARA), T(NPARA) DIMENSION COV(NPARA,NPARA), Y(2000,3), YNORM(2000) DIMENSION V(3,3), S(3,3), INDX(2000) DIMENSION B(2000) DIMENSION AA(2000,3), COVTR(3,3), ERROR(3,3) C---------------------------------------------------------------- C COMPUTES A ROBUST COVARIANCE MATRIX BASED ON C MULTIVARIATE TRIMMING (MVT) C C INPUT : A = DATA MATRIX C : IRDA = ROW DIMENSIONS IN CALLING PROGRAM C : NDATA= NUMBER OF DATA C : COV = STARTING COVARIANCE MATRIX (LOCAL) C : T = ESTIMATE OF LOCATION VECTOR C C OUTPUT : COVTR= ESTIMATE OF COVARIANCE MATRIX C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C----------------------------------------------------------------- ICOUNT = 0 EN = DFLOAT(NDATA) K = ALF * EN ITRIM = NDATA - K CALL COVMTX(A, IRDA, NDATA, NPARA,T, COVTR) 5 ICOUNT = ICOUNT + 1 CALL CHOLSK(COVTR, 3, 3, S, 3) CALL XNVLOW(S, 3, 3, V) CALL YOUT(A, IRDA, NDATA, NPARA, V, T, Y, YNORM) C---------------------------------------------------------------- C SORT DATA C---------------------------------------------------------------- NPL = 1 CALL SORTD(YNORM, NDATA, 1, NDATA, 1, INDX, NPL, B, # YNORM, YNORM, YNORM, YNORM, YNORM, YNORM, YNORM, YNORM,YNORM) DO 10 I =1, ITRIM J = INDX(I) AA(I,1) = A(J,1) AA(I,2) = A(J,2) AA(I,3) = A(J,3) 10 CONTINUE CALL COVMTX(AA, IRDA, ITRIM, NPARA,T, COV) C CONST = EN / ( DFLOAT(ITRIM) ) CONST = 1.D0 CALL CONMLT(COV, NPARA, CONST) C WRITE(6,200) ICOUNT 200 FORMAT(3X,'COVARIANCE MATRIX FROM MULTIVARIATE TRIMMING', # 3X,'ITERATION NUMBER=',I5) C--------------------------------------------------- C CHECK CONVERGENCE C--------------------------------------------------- CALL MSUBD(ERROR, 3, COV, 3, COVTR, 3, 3, 3) DE = FNORM(ERROR, 3) IF(DE .GT. 10.D-03) THEN DO 40 I=1,3 DO 50 J=1,3 COVTR(I,J) = COV(I,J) 50 CONTINUE 40 CONTINUE GO TO 5 ENDIF C CALL MOUTD(COV, 3, 3, 3) RETURN END C C SUBROUTINE GROSS(A, IRDA, NDATA, NPARA, AH, BH, CRIT, T, # V, COV, INDEX, ICLEAN) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(IRDA,NPARA), T(NPARA), V(NPARA,NPARA) DIMENSION COV(3,3) DIMENSION INDEX(NDATA) DIMENSION Y(2000,3), YNORM(2000) C-------------------------------------------------------------- C MULTIVARIATE OUTLIER DETECTION BASED ON CAMPBELL'S WEIGHTS C CREATES A NEW MATRIX (I.E.,CLEAN) WHICH CONTAINS C ACCEPTABLE OBSERVATIONS ONLY. C C INPUT : A = DATA MATRIX C : IRDA = ROW DIMENSIONS IN CALLING PROGRAM C : NDATA = NUMBER OF DATA POINTS C : NPARA = NUMBER OF PARAMETRES C : AH = CAMPBEEL'S CONSTANTS FOR R0 C : BH = CAMPBEEL'S CONSTANTS B2 C : CRIT = MAXIMUM WEIGHT ACCEPTABLE (CRITERION) C : V = LOWER TRIANGULAR MATRIX OF CHOLESKI DECOMP. C : COV = FINAL COVARIANCE MATRIX C C OUTPUT : INDEX = VECTOR CONTAINING INDICES OF CLEANED DATA C : A = MATRIX CONTAINING CLEANED DATA C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C-------------------------------------------------------------------- CALL YOUT(A, IRDA, NDATA, NPARA, V, T, Y, YNORM) J = 0 DO 10 I=1, NDATA R = YNORM(I) U = WSCAM(AH, BH, R) WI = U / R C WRITE(15 ,*) R, WI C WRITE(6 ,*) R, WI, I IF( WI .GT. CRIT) THEN J = J + 1 A(J,1) = A(I,1) A(J,2) = A(I,2) A(J,3) = A(I,3) INDEX(J) = I ENDIF 10 CONTINUE ICLEAN = J 100 FORMAT(2X,' DX=',F10.2,' DY=',F10.2,' DZ=',F10.2,' RI=',F10.2, # ' J=',I5,' I=',I5) RETURN END C FUNCTION R0CAM(NPARA, B1) IMPLICIT REAL * 8 (A-H,O-Z) C---------------------------------------------------------- C COMPUTES CAMPBELL'S WEIGTHING FUNCTION FOR THE MEAN C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C----------------------------------------------------------- R0CAM = DSQRT( DFLOAT(NPARA) ) + B1 / DSQRT(2.0D0) RETURN END FUNCTION WSCAM(AM0, B2, R) IMPLICIT REAL * 8 (A-H,O-Z) C---------------------------------------------------------- C COMPUTES CAMPBELL'S WEIGTHING FUNCTION FOR THE COVARIANCE C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C----------------------------------------------------------- ERROR = 0.1D-20 ALLOW = - DLOG( ERROR ) IF(R .LE. AM0) THEN WSCAM = R RETURN ENDIF IF(R .GT. AM0) THEN ARG =( (R - AM0) / B2 ) ** 2 TEST = ARG IF( TEST .GT. ALLOW) THEN TEMP = 0.0D0 GO TO 5 ENDIF ARG = - ARG / 2.0D0 TEMP = DEXP( ARG ) 5 WSCAM = AM0 * TEMP C WRITE(6,100) R, ARG, B2, AM0, WSCAM RETURN ENDIF 100 FORMAT(3X,'R=',F10.2,3X,'ARG=',F10.5,2X,'B2=',F10.4,' R0=', # F10.4,' OMEGA=',G10.4) END C FUNCTION CAMDRV(A, B, R) IMPLICIT REAL * 8 (A-H,O-Z) C---------------------------------------------------------- C COMPUTES DERIVATIVE OF C CAMPBELL'S WEIGTHING FUNCTION FOR THE COVARIANCE C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C----------------------------------------------------------- ERROR = 0.1D-20 ALLOW = - DLOG( ERROR ) IF(R .LE. A) THEN CAMDRV = 2.D0 * R RETURN ENDIF IF(R .GT. A) THEN ARG =( (R - A) / B ) ** 2 TEST = ARG IF( TEST .GT. ALLOW) THEN TEMP = 0.0D0 GO TO 5 ENDIF ARG = - ARG TEMP = DEXP( ARG ) 5 P = - (A / B ) ** 2 Q = 2.D0 * (R - A) CAMDRV = P * Q * TEMP C WRITE(6,100) R, ARG, B, A, CAMDRV RETURN ENDIF 100 FORMAT(3X,'R=',F10.2,3X,'ARG=',F10.5,2X,'B =',F10.4,' A =', # F10.4,' CAMDRV =',G10.4) END C SUBROUTINE LOCML(A, IRDA, NDATA, NPARA, V, AM, BM, # EPSI, T, H, VHNRM, MODE,IP) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(IRDA,NPARA), T(NPARA), H(NPARA) DIMENSION Y(2000,3), YNORM(2000), V(3,3), VH(3) C DO 90 K=1, 20 CALL YOUT(A, IRDA, NDATA, NPARA, V, T, Y, YNORM) CALL MLOCAT(A, IRDA, NDATA, NPARA, T, YNORM, AM, BM, H,MODE) DO 80 I= 1, NPARA T(I) = T(I) + H(I) 80 CONTINUE CALL MMULD(VH, IRDA, V, NPARA, H, NPARA, 3, 3, 1) CALL VNORM(VH, NPARA, VHNRM) IF(IP .GT. 0) THEN WRITE(6,110) WRITE(6,1000) CALL MOUTD(T,NPARA,3,1) WRITE(6,120) WRITE(6,1000) CALL MOUTD(H,NPARA,3,1) WRITE(6,100) VHNRM, K ENDIF IF(VHNRM .LT. EPSI) GO TO 91 90 CONTINUE C------------------------------------------------------------------- 100 FORMAT(2X,'NORM OF THE ERROR =',F10.7,' WITH ITERATION',I7) 110 FORMAT(10X,'LOCATION ESTIMATE T') 120 FORMAT(10X,' CORRECTION H ') 1000 FORMAT(5X,'------------------------------------------------------- #-----------------------------------') C-------------------------------------------------------------------- 91 RETURN END C SUBROUTINE COVML(A, IRDA, NDATA, NPARA, T, AH, BH, EPSI, V, COV #,IP,MODE) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(IRDA,NPARA), T(NPARA), WORK(2000) DIMENSION S(3,3), ST(3,3), B(3,3), C(3,3) DIMENSION SINV(3,3), IW1(3), IW2(3) DIMENSION V1(2000) DIMENSION ARRAY(2000), ARRAY1(2000) DIMENSION COV(3,3) , XI(3,1), DXI(3,1) C DIMENSION V(3,3), Y(2000,3), YNORM(2000) DIMENSION W(3,3) DIMENSION ALFA(3,1), BETA(1,3), OUT(3,3) DIMENSION COVUP(3,3), H(3), VH(3) C----------------------------------------------------------------- C COMPUTES AN M-ESTIMATE OF THE COVARIANCE MATRIX C COVARIANCE MAXIMUM LIKELIHOOD C A = DATA MATRIX (INPUT) C IRDA = ROW DIMENSIONS OF A MATRIX IN CALLING PROGRAM (INPUT) C NDATA= NUMBER OF DATA (INPUT) C NPARA= NUMBER OF PARAMETRES (INPUT) C T = ESTIMATE OF LOCATION VECTOR (INPUT) C AH = HUBER'S A CONSTANT USUALLY=0.0 (INPUT) C BH = HUBER'S B CONSTANT (INPUT) C EPSI = EPSILON ERROR USED FOR CONVERSION (INPUT) C V = LOWER TRIANGULAR MATRIX (VT * V =INV{COV}) (INPUT) C COV = FIMAL M-ESTIMATE OF COVARIANCE MARTIX (OUTPUT) C IP = 1 PRINTS OUT MATRICES INVOLVED IN COMPUTATIONS C = 0 NO PRINTING TAKES PLACE EXCEPT THE FINAL COVARIANCE C MODE = 1 HUBER'S APPROACH IS FOLLOWED C = 2 CAMPBELL'S APPROACH IS FOLLOWED C = 3 MARONNA'S APPROACH WITH CAUCHY DISTRIB. (D.F.=1) C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C------------------------------------------------------------------- C START WITH CLASSICAL ESTIMATES FOR COVARIANCE C-------------------------------------------------------------------- IF(IP .GT. 0) THEN WRITE(6,100) 100 FORMAT(6X,'..............INPUT COVARIANCE MATRIX....... # ................................',//) CALL MOUTD(COV, 3, 3, 3) ENDIF C------------------------------------------------------------------ C COMPUTE CHOLESKI DECOMPOSITION C------------------------------------------------------------------- CALL CHOLSK(COV, 3, 3, S,NPARA) IF(IP .GT. 0) THEN WRITE(6,110) 110 FORMAT(6X,'..............CHOLESKI DECOMPOSITION.............') CALL MOUTD(S, 3, 3, 3) ENDIF C---------------------------------------------------------------- C COMPUTE INVERSE OF LOWER TRIANGURAR C---------------------------------------------------------------- CALL XNVLOW(S, NPARA, NPARA, V) IF(IP .GT. 0) THEN WRITE(6,120) 120 FORMAT(6X,'............ V MATRIX ........................') CALL MOUTD(V, 3, 3, 3) ENDIF C---------------------------------------------------------------- C CONSTRUCT VECTORS Y=V(XI-T) C---------------------------------------------------------------- ICOUNT = 0 DO 1 I = 1,25 CALL YOUT(A, IRDA, NDATA, NPARA, V, T, Y, YNORM) C--------------------------------------------------------------- C COMPUTE M-ESTIMATE FOR THE COVARIANCE MATRIX C--------------------------------------------------------------- CALL MCOVAR(Y, IRDA, NDATA, NPARA, YNORM, T, AH, BH, # COVUP,MODE) IF(IP .GT. 0) THEN WRITE(6,130) 130 FORMAT(6X,'-----------C-MATRIX----------------------------') CALL MOUTD(COVUP, NPARA, NPARA, 3) ENDIF C--------------------------------------------------------------- C TAKE THE CHOLESKI DECOMPOSITION C--------------------------------------------------------------- CALL CHOLSK(COVUP, 3, 3, S,NPARA) IF(IP .GT. 0) THEN WRITE(6,140) 140 FORMAT(6X,'---------- CHOLESKI DECOMPOSITION-(M+1)--------') CALL MOUTD(S , NPARA, NPARA, 3) ENDIF C--------------------------------------------------------------- C COMPUTE THE INVERSE OF THE LOWER TRIANGULAR C--------------------------------------------------------------- CALL XNVLOW(S, NPARA, NPARA, W) IF(IP .GT. 0) THEN WRITE(6,150) 150 FORMAT(6X,'---------- W-MATRIX----------------------------') CALL MOUTD(W , NPARA, NPARA, 3) ENDIF C--------------------------------------------------------------- C COMPUTE FINAL ESTIMATE V:= W V C--------------------------------------------------------------- CALL MMULD(V, NPARA, W, NPARA, V, NPARA, 3, 3, 3) C--------------------------------------------------------------- C FINAL COVARIANCE M-ESTIMATE C--------------------------------------------------------------- CALL MTXMD(COV, NPARA, V, NPARA, 3, 3, 3) CALL MINVD(COV, NPARA, 3, DETA, IW1, IW2) IF(IP.GT.0)CALL MOUTD(COV , NPARA, NPARA, 3) C--------------------------------------------------------------- C COMPUTE THE NORM OF THE ERROR (W-I) C--------------------------------------------------------------- CALL ERROR(W,NPARA, WNORM) IF(IP .GT. 0) WRITE(6,170) WNORM, EPSI ICOUNT = ICOUNT + 1 IF(WNORM .LT. EPSI) GO TO 9999 1 CONTINUE 9999 IF(IP .GT. 0) THEN WRITE(6,200) AH,BH,ICOUNT,EPSI CALL MOUTD(COV , NPARA, NPARA, 3) ENDIF 200 FORMAT(2X,'FINAL M-ESTIMATE OF COVARIANCE MATRIX WITH AV=', #F5.2,2X,'BV=',F5.2,2X,'AFTER',I3,1X,'ITERATIONS',1X, #'WITH EPSILON=', #F10.7) 170 FORMAT(3X,'FROBENIUS NORM OF THE ERROR=',F10.7,2X,'EPSILON=', #F10.7) RETURN END C FUNCTION WBR(B, R) IMPLICIT REAL * 8(A-H,O-Z) C------COMPUTES HUBER'S WEIGHT FOR LOCATION-------- WBR = DMIN1(1.D0, B/R ) RETURN END C SUBROUTINE VNORM(V1, N, VNRM) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION V1( N ) SUM = 0.D0 DO 10 I=1,N SUM = SUM + V1(I) ** 2 10 CONTINUE VNRM = SUM RETURN END C SUBROUTINE MLOCAT(A, IRDA, NDATA, NPARA, T, YNORM, AM,BM, # H, MODE) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(IRDA,NPARA), T(NPARA), YNORM(NDATA), H(NPARA) DIMENSION DXI(3), ADDH(3) C---------------------------------------------------------------- C COMPUTES THE M-LOCATION ESTIMATE C C AUTHOR : S.MERTIKAS C DATE : NOVEMBER,1986 C--------------------------------------------------------------- ADDW = 0.0D0 ADDH(1) = 0.0D0 ADDH(2) = 0.0D0 ADDH(3) = 0.0D0 DO 10 I =1, NDATA R = YNORM(I) IF(MODE .EQ. 1) THEN WI = WBR(BM, R) ENDIF IF( MODE .EQ. 2) THEN U = WSCAM(AM, BM, R) WI = U / R ENDIF IF(MODE .EQ. 3) THEN WI = ( DFLOAT(NPARA) + 1.D0 ) / (1.D0 + R**2) ENDIF ADDW = ADDW + WI DO 20 J =1, NPARA DXI(J) = A(I,J) - T(J) H(J) = WI * DXI(J) ADDH(J)= ADDH(J) + H(J) 20 CONTINUE 10 CONTINUE DO 30 J = 1, NPARA H(J) = ADDH(J) / ADDW 30 CONTINUE RETURN END C SUBROUTINE ERROR(W, IRDA, WNORM) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION W(IRDA,IRDA), ONE(3,3) C------------------------------------------------- C COMPUTES THE FROBENIUS NORM OF THE (W-I) MATRIX C------------------------------------------------- SUM = 0.0D0 DO 10 I=1, IRDA DO 11 J = 1,IRDA ONE(I,J) = 0.0D0 11 CONTINUE 10 CONTINUE DO 12 I=1,IRDA ONE(I,I) = -1.D0 12 CONTINUE CALL MADDD(ONE, IRDA, W, IRDA, ONE, IRDA, 3, 3, 3) DO 20 I =1, IRDA DO 30 J =1, IRDA SUM = SUM + ONE(I,J)**2 30 CONTINUE 20 CONTINUE WNORM = DSQRT( SUM ) 999 RETURN END C FUNCTION FNORM(W, N) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION W(N,N) C------------------------------------------------------- C COMPUTES THE FORBENIUS NORM OF A MATRIX C------------------------------------------------------- SUM = 0.0D0 DO 10 I =1, N DO 20 J =1, N SUM = SUM + W(I,J) ** 2 20 CONTINUE 10 CONTINUE FNORM = DSQRT( SUM ) RETURN END C SUBROUTINE MCOVAR(Y, IRDA, NDATA, NPARA, YNORM, T, AH, BH, # COVUP,MODE) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION Y(IRDA, NPARA), YNORM(NDATA), T(NPARA) DIMENSION COVUP(NPARA, NPARA) DIMENSION ALFA(3,1), BETA(1,3), OUT(3,3) C--------------------------------------------------------------------- C COMPUTES THE M-ESTIMATE FOR THE COVARIANCE MATRIX USING C HUBER'S AND MARONNA'S ALGORITHMS C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C--------------------------------------------------------------------- ADDS = 0.D0 DO 1 I=1,NPARA DO 2 J=1,NPARA COVUP(I,J) = 0.D0 2 CONTINUE 1 CONTINUE C DO 10 I= 1, NDATA DO 20 J=1, NPARA ALFA(J,1) = Y(I,J) 20 CONTINUE CALL OUTPRD(ALFA, NPARA, OUT, NPARA) IF( MODE .EQ. 1) THEN R = YNORM(I) U = UBR(AH, BH, R) S = SHUBER(U, R) ENDIF C WRITE(6,500) R,U,S , I 500 FORMAT(2X,'R=',F10.2,3X,'U=',F12.4,3X,'S=',F12.4,3X,'I=',I5) IF(MODE .EQ. 2) THEN R = YNORM(I) U = WSCAM(AH, BH, R) S = ( U / R ) ** 2 C WRITE(6,100) R, U, S, I 100 FORMAT(3X,'R=',F10.2,' OMEGA=',F10.4,' S=',F12.6,' I=',I6) ENDIF IF(MODE .EQ. 3) THEN R = YNORM(I) S = (DFLOAT(NPARA) + 1.D0) / (1.D0 + R**2) ENDIF CALL CONMLT(OUT, NPARA, S) ADDS = ADDS + S CALL MADDD(COVUP,NPARA, OUT, NPARA, COVUP, NPARA, NPARA, NPARA) 10 CONTINUE IF( MODE .EQ. 2 ) ADDS = ADDS -1.D0 IF( MODE .EQ. 1 .OR. MODE .EQ. 3 ) ADDS = DFLOAT(NDATA) CONST = 1.D0 / ADDS CALL CONMLT(COVUP, NPARA, CONST) RETURN END C SUBROUTINE YOUT(A, IRDA, NDATA, NPARA, V, T, Y, YNORM) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION A(IRDA, NPARA), V(NPARA, NPARA), T(NPARA) DIMENSION Y(IRDA, NPARA), YNORM(NDATA), YI(3,1), DXI(3,1) C-------------------------------------------------------------------- C COMPUTES THE Y VECTOR (=V(XI-T)) AND THE NORM OF YI C C A = DATA MATRIX (INPUT) C IRDA = ROW DIMENSION IN CALLING PROGRAM (INPUT) C NDATA= NUMBER OF DATA (INPUT) C NPARA= NUMBER OF PARAMETRES (INPUT) C V = INVERSE OF THE LOWER TRIANGULAR MATRIX (INPUT) C T = ESTIMATE OF LOCATION VECTOR (INPUT) C Y = MATRIX CONTAINING ALL YI'S (OUTPUT) C YNORM= VECTOR CONTAINING ALL SQUARE-ROOT NORMS OF YI'S (OUTPUT) C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C--------------------------------------------------------------------- DO 9 I=1,NDATA YNORM(I) = 0.D0 9 CONTINUE DO 10 I = 1, NDATA DO 20 J =1, NPARA DXI(J,1) = A(I,J) - T(J) 20 CONTINUE CALL MMULD(YI, NPARA, V, NPARA, DXI, NPARA, 3, 3, 1) CALL VNORM(YI, NPARA, R2) YNORM(I) = DSQRT( R2 ) Y(I,1) = YI(1,1) Y(I,2) = YI(2,1) Y(I,3) = YI(3,1) 10 CONTINUE RETURN END C FUNCTION UBR(AH, BH, R) IMPLICIT REAL * 8 (A-H,O-Z) C--------- A AND B ARE HUBER'S CONSTANTS AND R IS THE DISTANCE IF(R .LE. AH) THEN UBR = AH**2 RETURN ENDIF IF(R .GE. AH .AND. R .LE. BH) THEN UBR= R**2 RETURN ENDIF UBR= BH**2 RETURN END C FUNCTION SHUBER(U, R) IMPLICIT REAL * 8 (A-H,O-Z) C--------- COMPUTES THE S-WEIGHT FUNCTION--------------------- SHUBER = U / (R**2) RETURN END C SUBROUTINE OUTPRD(ALFA, N, OUT, IRDC) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION ALFA(IRDC,1), BETA(1,3), OUT (IRDC,N) C----------COMPUTES OUTER PRODUCT OF ALFA AND BETA VECTORS------- CALL TRNSD(BETA, 1, ALFA, N, N, 1) CALL MMULD(OUT, IRDC, ALFA, IRDC, BETA, 1, N, 1, N) RETURN END C SUBROUTINE COVMTX(A, IRDA, NDATA, NPARA,T, COV) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION A(IRDA,NPARA), T(NPARA), COV(NPARA, NPARA) DIMENSION XI(3,1), DXI(3,1), DXITR(1,3), ADD(3,3) C-------------------------------------------------------------- C COMPUTES AN ESTIMATE OF A COVARIANCE MATRIX C C A = DATA MATRIX (INPUT) C IRDA = ROW DIMENSIONS OF A IN CALLING PROGRAM (INPUT) C NDATA = NUMBER OF DATA (INPUT) C NPARA = NUMBER OF PARAMETRES TO BE ESTIMATED (INPUT) C T = ESTIMATE OF LOCATION VECTOR (INPUT) C COV = ESTIMATE OF COVARIANCE MATRIX (OUTPUT) C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C-------------------------------------------------------------- C C NVEC = 3 IDC = NPARA C DO 6 I =1, NPARA DO 5 J =1, NPARA COV(I,J) = 0.D0 ADD(I,J) = 0.D0 5 CONTINUE 6 CONTINUE C DO 10 I=1, NDATA DO 20 J =1, NPARA XI(J,1) = A(I,J) 20 CONTINUE C CALL SUBTRA(XI, T, NPARA, DXI) CALL TRNSD(DXITR, 1, DXI, NVEC, NPARA, 1) CALL MMULD(ADD, IDC, DXI, NVEC, DXITR, 1, NPARA, 1, NPARA) CALL MADDD(COV, IDC, COV, IDC, ADD, IDC, NPARA, NPARA) 10 CONTINUE CONST = 1.D0 / DFLOAT( NDATA-1 ) CALL CONMLT(COV, NPARA, CONST) RETURN END C SUBROUTINE CONMLT(COV, NPARA, CONST) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION COV(NPARA, NPARA) C--------------------------------------------------------- C MULTIPLIES A COVARIANCE MATRIX BY A CONSTANT C--------------------------------------------------------- DO 10 I =1, NPARA DO 20 J=1, NPARA COV(I,J) = CONST * COV(I,J) 20 CONTINUE 10 CONTINUE RETURN END C SUBROUTINE SUBTRA(ARRAY, T, NPARA, ARRAY1) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION ARRAY(NPARA), T(NPARA), ARRAY1(NPARA) C--------------------------------------------------------- C SUBTRACTS TWO VECTORS (THAT IS, T FROM XI) C--------------------------------------------------------- DO 10 I = 1, NPARA ARRAY1(I) = ARRAY(I) - T(I) 10 CONTINUE RETURN END C SUBROUTINE XNVLOW(S, IRDS, N, SINV) IMPLICIT REAL * 8(A-H,O-Z) DIMENSION S(IRDS, N), SINV(IRDS,N) C------------------------------------------------------------------ C COMPUTES THE INVERSE OF A LOWER TRIANGULAR MATRIX S C C S = LOWER TRIANGULAR MATRIX TO BE INVERTED C IRDS = DECLARED DIMENSIONS IN CALLING PROGRAM C N = DIMENSIONS IN THIS PROGRAM C SINV = INVERTED MATRIX (OUTPUT) C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C----------------------------------------------------------------- DO 10 I=1, N SINV(I,I) = 1.D0 / S(I,I) 10 CONTINUE SINV(1,2) = 0.D0 SINV(1,3) = 0.D0 SINV(2,3) = 0.D0 C DO 20 J = 1, N - 1 DO 30 I = J+1, N SUM = 0.0D0 DO 40 K= J, I-1 SUM = SUM + S(I,K) * SINV(K,J) 40 CONTINUE SINV(I,J) = - SUM / S(I,I) 30 CONTINUE 20 CONTINUE RETURN END C SUBROUTINE CHOLSK(A, IRDA, NA, S, IDC) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(IRDA,NA), S(IDC,NA) C---------------------------------------------------------------------- C COMPUTES THE CHOLESKI DECOMPOSITION OF MATRIX A( = S * TRANS(S)) C A = MATRIX (IRDA,NA) POSITIVE DEFINITE AND SYMMETRIC C IRDA = ROW DIMENSION OF A C NA = RANK OF MATRIX A C S = OUTPUT MATRIX (LOWER TRIAGULAR) C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 C--------------------------------------------------------------------- IF(NA .LT. 1) GO TO 999 CALL CONMLT(S, 3, 0.D0) S(1,1) = DSQRT( A(1,1) ) DO 10 I = 2,NA S(I,1) = A(I,1) / S(1,1) 10 CONTINUE DO 50 J = 2, NA SUM = 0.0D0 J1 = J - 1 DO 20 K = 1, J1 SUM = SUM + S(J,K) ** 2 20 CONTINUE S(J,J) = DSQRT( A(J,J) - SUM ) IF( J .EQ. NA ) GO TO 50 J2 = J + 1 DO 40 I = J2, NA SUM = 0.D0 DO 30 K = 1, J1 30 SUM = SUM + S(I,K)* S(J,K) S(I,J) = ( A(I,J) - SUM) / S(J,J) 40 CONTINUE 50 CONTINUE RETURN 999 WRITE(6,1000) 1000 FORMAT(4X,'DIMENSIONS OF A MATRIX ARE LESS THAN 1 --SUBROUTINE # CHOLSK',//) END C SUBROUTINE LOCATE(A, IRDA, NDATA, NPARA, T) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(IRDA, NPARA), T(NPARA), WORK(2000) C------------------------------------------------------------------ C COMPUTES THE LOCATION VECTOR BASED ON MEDIAN C A = DATA MATRIX(IN) C IRDA = NUMBER OF ROWS IN CALLING PROGRAM C NPARA = NUMBER OF PARAMETRES (IN) C NDATA = NUMBER OF DATA (IN) C T = ESTIMATE OF LOCATION VECTOR USING MEDIAN C C AUTHOR : S. MERTIKAS C DATE : NOVEMBER,1986 (FIRST SNOW IN FREDERICTON TODAY !!!!!) C------------------------------------------------------------------- DO 10 I=1, NPARA IWISH = I CALL SET(A, IRDA, NDATA, NPARA, IWISH, WORK) CALL SORT8(WORK, NDATA, IERROR) C N = NDATA C CALL XMEDIN(WORK, NDATA,XMEDN) CALL XMEAN(WORK, NDATA, XMEDN) T(I) = XMEDN 10 CONTINUE RETURN END C SUBROUTINE SET(A, IRDA, NDATA, NPARA, IWISH, WORK) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(IRDA, NPARA), WORK(NDATA) C DO 10 I=1,NDATA WORK(I) = A(I, IWISH) 10 CONTINUE RETURN END C SUBROUTINE SORT8(V1,N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION V1(N) I = 1 1 I = I + I IF(I.LE.N) GO TO 1 M = I - 1 2 M = M / 2 IF(M.EQ.0) RETURN K = N - M DO 4 J=1,K L = J 5 IF(L .LT. 1) GO TO 4 IF(V1(L+M) .GE. V1(L)) GO TO 4 X = V1(L+M) V1(L+M) = V1(L) V1(L) = X L = L - M GO TO 5 4 CONTINUE GO TO 2 END C SUBROUTINE XMEDIN(ARRAY,N,XMEDN) IMPLICIT REAL*8(A-H,O-Z) DIMENSION ARRAY(N) M = (N +1) / 2 L = (N +2) / 2 XMEDN = .5D0*(ARRAY(M) + ARRAY(L)) RETURN END //* //LKED.USERLIB DD DSN=A.M12129.SELIBOBJ,DISP=SHR //GO.FT15F001 DD DSN=&&TEMPPVA,DISP=(NEW,PASS), // UNIT=DASD,SPACE=(TRK,(1,1),RLSE), // DCB=(RECFM=FB,LRECL=80,BLKSIZE=8000) //GO.FT22F001 DD DSN=A.M12128.MERTIKAS.FINAL, // UNIT=P3350,DISP=OLD //*GO.FT17F001 DD DSN=A.M12128.MERTIKAS.GEDOP, //* UNIT=P3350,DISP=OLD //GO.SYSIN DD * -7.75 5.81 20 14.04 25.28 20 8.28 20.13 20 -20.05 3.83 20 20.28 40.03 20 23.29 42.53 20 35.99 -22.62 20 -4.26 0.51 20 10.77 13.93 20 -1.49 2.79 20 -7.34 -0.94 20 13.97 17.90 20 -2.15 0.52 20 2.36 4.55 20 2.60 4.95 20 7.69 12.49 21 //