C*********************************************************************** C SUBROUTINE PEREST HANDLES ALL ESTIMATION FROM PERPENDICULAR DIST- C ANCES FOR BOTH GROUPED AND UNGROUPED DATA. IT CALLS THE PROPER C ESTIMATION ROUTINES FOR NPD ESTIMATORS TO CALCULATE F(0) AND ITS C VARIANCE. IT THEN CALLS PARAM TO CALCULATE A DENSITY ESTIMATE AND C TO OUTPUT A VARIETY OF STATISTICS. THEN IT PLOTS THE FIT OF THE C ESTIMATED FUNCTION TO THE HISTOGRAM AND CALLS TEST TO DO PAR CHI- C SQUARE GOODNESS OF FIT TEST. FOR THE NEGATIVE EXPONENTIAL ESTI- C MATOR IT ALSO CALLS EXPT. C C SUBROUTINES CALLED: HEADER,MLE,FSEST,PARAM,NAME,EXPT,SETUP,PLT C PRNT,CLEAR,DVSION,PLTMOD,CELLS,TEST C*********************************************************************** SUBROUTINE PEREST C*********************************************************************** C DECLARATIONS C*********************************************************************** INCLUDE 'PARMTR.INC' LOGICAL GRP, POOL, PEST, SEST, DESC, DEF, TRUNC, CUTP, WARN, TR, 1 ANY LOGICAL SET, STRT, HELP, SK INTEGER PDEST, SDEST, CNT, SYSIN, FILPOS CHARACTER*1 HEAD, LABEL, UL, XT1(36), YT1(36), CHAR, CHARF CHARACTER*66 FILEIN,FILOUT,FILDOC INTEGER NO(MAXCEL) DOUBLE PRECISION PCELL, CELL REAL LOGL, MAX, FQ(MAXCEL), MESANG, EFREQ(MAXCEL) C*********************************************************************** C COMMON STATEMENTS C*********************************************************************** COMMON /TS/ Z(MAXOBJ) COMMON /PDOPT/ STRT, SET, PSTRT(10,MAXPAR), NPSET(10) COMMON /ALPHA/ LABEL(80), HEAD(30), UL(3,25) COMMON /PAGE/ IPAGE COMMON /INTER/ KCUT, CUT(MAXCEL), FREQ(MAXCEL), NCUT, NKC(5), 1 NK(5), RFREQ(5,MAXCEL), CCUT(5,MAXCEL) COMMON /ESTM/ PDEST(10), SDEST(3), NSD, NPD COMMON /OPTION/ GRP, POOL, PEST, SEST, DESC, DEF, CUTP, TRUNC, 1 HELP COMMON /STORE/ D(13), DL(13), DU(13), STD(13), LOGL(10), PB(5,10) COMMON /NUM/ XL(MAXLIN), WIDTH, N, CNT, CONV(3), VARN, IDF, WARN COMMON /MEASUR/ PD(MAXOBJ), SD(MAXOBJ), MESANG(MAXOBJ), 1 COMANG(MAXOBJ), SINA(MAXOBJ) COMMON /SOLN/ FZERO, VARF, FMAX, FMIN, FZ(100) COMMON /PART/ PCELL(MAXCEL,MAXPAR), CELL(MAXCEL) DOUBLE PRECISION PAR, VCMAT, G, XLL COMMON /DPAR/ PAR(MAXPR2), VCMAT(MAXPAR,MAXPAR), G(MAXPAR), 1 XLL, NPAR, INDEX COMMON /FILE/ SYSIN, SK, FILPOS COMMON /FILES/ FILEIN,FILOUT,FILDOC INCLUDE 'SCREEN.INC' C*********************************************************************** C DATA STATEMENTS C*********************************************************************** DATA CHAR /'*'/, CHARF /'f'/ DATA YT1 /9*' ','P','r','o','b','a','b','i','l','i','t','y',' ', 1 'D','e','n','s','i','t','y',8*' '/ DATA XT1 /6*' ','P','e','r','p','e','n','d','i','c','u','l','a', 1 'r',' ','D','i','s','t','a','n','c','e',8*' '/ C*********************************************************************** C PRINT OUT HEADING AND EXPLANATORY INFORMATION C*********************************************************************** CALL HEADER (1) WRITE (6,300) CALL HEADER (2) IF (HELP) THEN WRITE (6,510) FILDOC(FILPOS:)='NARRAT2.DOC' OPEN(UNIT=1,FILE=FILDOC,STATUS='OLD',ERR=20) GO TO 25 20 WRITE(0,*) 'Cannot open documentation file ',FILDOC STOP 'Documentation file error' 25 DO 30 I=1,49 READ(UNIT=1,FMT='(A79)') LINE 30 WRITE(UNIT=6,FMT='(1X,A79)') LINE CLOSE (UNIT=1) CALL HEADER (1) ENDIF WRITE (6,480) ANY=.FALSE. DO 40 I=1,NPD INDEX=PDEST(I) IF (INDEX.EQ.1) THEN WRITE (6,310) IF ((.NOT.SET).OR.(NPSET(I).EQ.0)) GO TO 40 WRITE (6,460) NPSET(I) ELSE IF (INDEX.EQ.2) WRITE (6,320) IF (INDEX.EQ.3) WRITE (6,330) IF (INDEX.EQ.4) WRITE (6,340) IF (INDEX.EQ.5) WRITE (6,360) NPSETI=NPSET(I) ENDIF IF (STRT.AND.(NPSET(I).NE.0)) WRITE (6,470) (PSTRT(I,J),J=1,NPSETI 1 ) IF ((INDEX.EQ.2).OR.(INDEX.EQ.3).OR.(INDEX.EQ.4).OR.(INDEX.EQ.5)) 1 ANY=.TRUE. 40 CONTINUE IF (ANY) WRITE (6,380) IF (WARN) VARN=FLOAT(CNT) XN=FLOAT(CNT) DO 290 I=1,NPD INDEX=PDEST(I) C*********************************************************************** C CALL ESTIMATION SUBROUTINE FOR GROUPED DATA C*********************************************************************** IF (GRP) THEN DO 50 J=1,KCUT 50 FREQ(J)=RFREQ(1,J) CALL MLE (I) ELSE C*********************************************************************** C CALL ESTIMATION SUBROUTINE FOR UNGROUPED DATA C*********************************************************************** IF (INDEX.EQ.1) CALL FSEST (I) IF (INDEX.EQ.2 .OR. INDEX.EQ.3) CALL MLE (I) IF (INDEX.NE.4) GO TO 90 C*********************************************************************** C FOR THE NEGATIVE EXPONENTIAL ESTIMATOR WITH UNTRUNCATED DATA PAR C CLOSED FORM SOLUTION EXISTS, SO IT IS CALCULATED HERE. C*********************************************************************** IF (TRUNC) GO TO 80 SUM=0.0 DO 70 J=1,CNT 70 SUM=SUM+PD(J) PAR(1)=XN/SUM FZERO=PAR(1) VCMAT(1,1)=PAR(1)*PAR(1)/XN VARF=VCMAT(1,1) NPAR=1 XLL=XN*DLOG(PAR(1))-PAR(1)*SUM LOGL(I)=XLL CALL HEADER (1) CALL NAME (INDEX) WRITE (6,490) XLL CALL HEADER (1) CALL NAME (INDEX) GO TO 120 C*********************************************************************** C FOR THE NEGATIVE EXPONENTIAL ESTIMATOR WITH TRUNCATED DATA AN C ITERATIVE SOLUTION MUST BE CALCULATED. C*********************************************************************** 80 CALL MLE (I) GO TO 120 90 IF (INDEX.NE.5) GO TO 120 C*********************************************************************** C FOR THE HALF-NORMAL ESTIMATOR FOR UNTRUNCATED DATA PAR CLOSED FORM C SOLUTION EXISTS AND IS CALCULATED HERE. C*********************************************************************** IF (TRUNC) GO TO 110 SUM=0.0 DO 100 J=1,CNT 100 SUM=SUM+PD(J)*PD(J) PAR(1)=XN/(2.*SUM) FZERO=2.*DSQRT(PAR(1))/SQRT(PI) VARF=2.*PAR(1)/(XN*PI) VCMAT(1,1)=2.*PAR(1)*PAR(1)/XN XLL=-PAR(1)*SUM+XN*ALOG(FZERO) NPAR=1 LOGL(I)=XLL CALL HEADER (1) CALL NAME (INDEX) WRITE (6,500) XLL CALL HEADER (1) CALL NAME (INDEX) GO TO 120 C*********************************************************************** C FOR THE HALF-NORMAL ESTIMATOR FOR TRUNCATED DATA REQUIRES AN C ITERATIVE SOLUTION TO BE CALCULATED. C*********************************************************************** 110 CALL MLE (I) ENDIF C*********************************************************************** C CALCULATE DENSITY AND ASSOCIATED STATISTICS AND WRITE THEM OUT C*********************************************************************** 120 CALL HEADER (2) NPSET(I)=NPAR IF (CNT.LT.N) WRITE (6,450) CNT IF ((INDEX.EQ.1.OR.INDEX.EQ.2).AND.(.NOT.TRUNC)) WRITE (6,440) 1 WIDTH CALL PARAM (XN,DNSTY,VARD,DLL,DUL) D(I)=DNSTY STD(I)=VARD DL(I)=DLL DU(I)=DUL C*********************************************************************** C FOR UNTRUNCATED DATA IF THE NEGATIVE EXPONENTIAL ESTIMATOR IS C REQUESTED THEN TRANSFORM THE DATA TO UNIFORM RANDOM VARIABLES AND C TEST. C*********************************************************************** IF (TRUNC.OR.GRP.OR.(INDEX.NE.4)) GO TO 140 DO 130 JJ=1,N 130 Z(JJ)=PD(JJ) CALL HEADER (1) CALL EXPT (Z,N) C*********************************************************************** C FOR EACH SET OF CUT POINTS PRODUCE PAR PLOT OF THE ESTIMATED C FUNCTION OVERLAYING THE HISTOGRAM SHAPE AND PAR CHI-SQUARE GOODNESS C OF FIT TEST FOR GROUPED DATA THERE IS ONLY ONE SET OF CUT POINTS. C*********************************************************************** 140 DO 280 KK=1,NCUT CALL SETUP CALL CLEAR IF (GRP) GO TO 160 C*********************************************************************** C RETRIEVE CUT POINTS AND NUMBERS OF OBSERVATIONS FOR THE KKTH C REPLICATE. C*********************************************************************** KCUT=NKC(KK) DO 150 J=1,KCUT CUT(J)=CCUT(KK,J) 150 FREQ(J)=RFREQ(KK,J) C*********************************************************************** C CALCULATE FREQUENCIES AND WIDTHS OF EACH INTERVAL C*********************************************************************** 160 HOLD=0.0 CALL PLTMOD MAX=0.0 DO 170 J=1,KCUT FQ(J)=FREQ(J)/(FLOAT(CNT)*(CUT(J)-HOLD)) NO(J)=IFIX((CUT(J)*65./CUT(KCUT))+.5) IF (MAX.LT.FQ(J)) MAX=FQ(J) HOLD=CUT(J) 170 CONTINUE C*********************************************************************** C FIND SCALE FOR PLOT C*********************************************************************** IF (MAX.LT.FMAX) MAX=FMAX MAX=1.1*MAX CALL DVSION (0.,CUT(KCUT),FMIN,MAX) XINC=CUT(KCUT)/65. YINC=(MAX-FMIN)/45. C*********************************************************************** C PLOT TOPS OF THE HISTOGRAM C*********************************************************************** IHOLD=1 DO 190 J=1,KCUT Y=FLOAT(IFIX(FQ(J)/YINC+.5))*YINC NOJ=NO(J) DO 180 K=IHOLD,NOJ X=FLOAT(K)*CUT(KCUT)/65. 180 CALL PLT (CHAR,X,Y) 190 IHOLD=NO(J)+1 C*********************************************************************** C PLOT SIDES OF HISTOGRAM C*********************************************************************** DO 220 J=1,KCUT NO2=IFIX((FQ(J)/YINC)+.5) IF (J.EQ.KCUT) GO TO 200 JP1=J+1 NO3=IFIX((FQ(JP1)/YINC)+.5) NO2=MAX0(NO2,NO3) 200 IF (NO2.EQ.0) GO TO 220 X=FLOAT(NO(J))*XINC DO 210 K=1,NO2 Y=FLOAT(K)*YINC 210 CALL PLT (CHAR,X,Y) 220 CONTINUE C*********************************************************************** C PLOT ESTIMATED FUNCTION C*********************************************************************** DO 230 J=1,65 X=FLOAT(J)*XINC 230 CALL PLT (CHARF,X,FZ(J)) C*********************************************************************** C PRINT THE PLOT C*********************************************************************** CALL HEADER (1) WRITE (6,350) KK CALL PRNT (XT1,YT1) WRITE (6,430) (UL(1,J),J=1,25) C*********************************************************************** C CHI SQUARE GOODNESS OF FIT TEST C*********************************************************************** CALL HEADER (1) WRITE (6,370) KK C*********************************************************************** C FOR THE FOURIER SERIES AND THE EXPONENTIAL POWER SERIES THE DATA C ARE CONSIDERED TRUNCATED BECAUSE PAR WIDTH MUST BE CHOSEN FOR EST- C IMATION. IF THE DATA ARE SPECIFIED UNTRUNCATED THEN THE LARGEST C OBSERVATION IS USED FOR THE WIDTH. C*********************************************************************** TR=.FALSE. IF (TRUNC.OR.(INDEX.EQ.1).OR.(INDEX.EQ.2)) TR=.TRUE. IF (GRP) GO TO 260 C*********************************************************************** C FOR UNGROUPED DATA THE FINAL CUT POINT MAY BE LESS THAN THE WIDTH C FOR HISTOGRAMS AND PLOTS BUT FOR CHI-SQUARE TESTS ALL OF THE DATA C USED FOR ESTIMATION MUST BE USED FOR TESTING. THIS IS ADJUSTED C HERE. C*********************************************************************** IF ((CUT(KCUT).EQ.WIDTH).OR.(.NOT.CUTP)) GO TO 250 SUM=0.0 DO 240 IJ=1,KCUT 240 SUM=SUM+FREQ(IJ) KCUT=KCUT+1 CUT(KCUT)=WIDTH FREQ(KCUT)=FLOAT(CNT)-SUM C*********************************************************************** C CALL THE ROUTINE WHICH CALCULATES THE EXPECTED FREQUENCIES FOR THE C FIT. C*********************************************************************** 250 CALL CELLS (2) 260 DO 270 JJ=1,KCUT 270 EFREQ(JJ)=SNGL(CELL(JJ)) CALL TEST (FREQ,EFREQ,KCUT,NPAR,XN,CUT,TR,P) PB(KK,I)=P 280 CONTINUE 290 CONTINUE C*********************************************************************** C FORMAT STATEMENTS C*********************************************************************** RETURN 300 FORMAT (//10X,60('*')/10X,'*',58X,'*'/10X,'*',6X, 1 'Density Estimation from Perpendicular Distances',5X,'*'/ 2 10X,'*',58X,'*'/10X,60('*')) 310 FORMAT (/6X,'Fourier Series (FSER)'/6X, 1'F(X) = 1/W + A(1)*COS(3.14159*X/W) + A(2)*COS(2*3.14159*X/W)'/ 213X,'+ ... + A(NPAR)*COS(NPAR*3.14159*X/W)'/ 313X,'NPAR is the number of parameters used in the estimator') 320 FORMAT (/6X,'Exponential power series (EXPS)'/ 1 6X,'G(X) = EXP(-(X/PAR(1))**PAR(2)), F(X) = G(X)/NORM') 330 FORMAT (/6X,'Exponential Polynomial (EXPL)'/ 1 6X,'G(X) = EXP(-PAR(1)*X - PAR(2)*X**2), F(X) = G(X)/NORM') 340 FORMAT (/6X,'Negative Exponential (NEXP)'/ 1 6X,'G(X) = EXP(-PAR(1)*X), F(X) = G(X)/NORM') 350 FORMAT ('0Estimated PDF of perpendicular distances (f) fit to', 1' histogram of data (*)'/ 2' Cut Point Set ',I1) 360 FORMAT (/6X,'Half-Normal (HNOR)'/ 1 6X,'G(X) = EXP(-PAR(1)*X**2), F(X) = G(X)/NORM') 370 FORMAT ('0Chi-square goodness of fit test of the null hypothesis', 1' that the'/ 2' model provides an adequate fit to the perpendicular', 2' distance data.'/'0Cut Point Set ',I1) 380 FORMAT ('0NORM represents the normalizing constant equal to the'/ 1' integral from 0 to W of G(X)') 430 FORMAT (35X,'in ',25A1) 440 FORMAT (' ***NOTE - The data are untruncated but a width must', 1' be chosen for the estimator.'/ 2' The last order statistic (largest measurement) was used'/ 3' for the WIDTH = ',G10.4) 450 FORMAT (' *** NOTE - The data have been truncated at a width'/ 113X,'such that some of the data have not been used.'/ 213X,'Only ',I4,' data points have been used.') 460 FORMAT (6X,'Number of parameters requested, NPAR = ',I2) 470 FORMAT (6X,'Starting values for maximum likelihood iteration:'/ 1 6X,8G9.3) 480 FORMAT (//'0 The following is a list of the estimators', 1' requested and representa-'/ 2' tions of their functional forms. If the number', 3' of parameters has been'/ 4' set or starting values have been given, then these', 4' are also listed.'//) 490 FORMAT (//'0The untruncated version of the negative exponential', 1' estimator does'/ 2' not require iterative techniques to find the solution. It can', 3' be solved in'/ 4' closed form. The log-likelihood value for this model (for use', 5' in comparison'/ 6' to other models is ',G11.5,' The estimate of density given is', 7' the true'/ 8' maximum likelihood estimate which is biased by a factor of', 9' (n-1)/n. This'/ A' is not a substantial bias for reasonable sample sizes. If', B' desired the'/ C' unbiased answer can be computed easily.') 500 FORMAT (//'0The untruncated version of the half-normal estimator', 1' does not require'/ 2' iterative techniques to find the solution. It can be solved', 3' in closed'/ 4' form. The log-likelihood value for this model (for use in', 5' comparison to'/ 6' other models is ',G11.5,' The estimate of density given is', 7' the true'/ 8' maximum likelihood estimate which is biased. However, at', 9' reasonable sample'/ A' sizes the bias is negligible. Quinn(1977) has given an', B' unbiased estimator.') 510 FORMAT (///) END