C*********************************************************************** C THIS SUBROUTINE HANDLES ALL ESTIMATION FOR UNGROUPED SIGHTING C DISTANCES AND ANGLES. IT CALLS THE PROPER ESTIMATION ROUTINE C OR DIRECTLY CALCULATES C(0) AND THEN CALCULATES F(0). IT CALLS C PARAM TO CALCULATE DENSITY AND ITS ASSOCIATED STATISTICS. C C SUBROUTINES CALLED: HEADER,ELLIPS,PARAM,SORT,HIST,CHI,UTEST C********************************************************************** SUBROUTINE SIGEST C*********************************************************************** C DECLARATIONS C*********************************************************************** INCLUDE 'PARMTR.INC' INTEGER SDEST, PDEST, CNT, STATUS, SYSIN, FILPOS LOGICAL WARN, STRT, MSET, GRP, POOL, PEST, SEST, DESC, DEF, CUTP, 1 TRUNC, SK LOGICAL HELP REAL SCUT(15), FR(15), MESANG, LOGL CHARACTER*66 FILEIN,FILOUT,FILDOC C*********************************************************************** C COMMON STATEMENTS C*********************************************************************** DOUBLE PRECISION PAR, VCMAT, G, XLL COMMON /DPAR/ PAR(MAXPR2), VCMAT(MAXPAR,MAXPAR), G(MAXPAR), 1 XLL, NPAR, INDEX COMMON /NUM/ XL(MAXLIN), WIDTH, N, CNT, CONV(3), VARN, IDF, WARN COMMON /IND/ IC, II, IREP, STATUS COMMON /PDOPT/ STRT, MSET, PSTRT(10,MAXPAR), NPSET(10) COMMON /ESTM/ PDEST(10), SDEST(3), NSD, NPD COMMON /SOLN/ FZERO, VARF, FMAX, FMIN, FZ(100) COMMON /MEASUR/ PD(MAXOBJ), SD(MAXOBJ), MESANG(MAXOBJ), 1 COMANG(MAXOBJ), SINA(MAXOBJ) COMMON /TS/ U(MAXOBJ) COMMON /STORE/ D(13), DL(13), DU(13), STD(13), LOGL(10), PB(5,10) COMMON /OPTION/ GRP, POOL, PEST, SEST, DESC, DEF, CUTP, TRUNC, 1 HELP COMMON /FILE/ SYSIN, SK, FILPOS COMMON /FILES/ FILEIN,FILOUT,FILDOC INCLUDE 'SCREEN.INC' IF (HELP) THEN FILDOC(FILPOS:)='NARRAT3.DOC' OPEN(UNIT=1,FILE=FILDOC,STATUS='OLD',ERR=3) GO TO 4 3 WRITE(0,*) 'Cannot open documentation file ',FILDOC GO TO 10 4 CALL HEADER (1) WRITE (6,160) WRITE(6,'(//)') DO 5 I=1,40 READ(UNIT=1,FMT='(A79)') LINE 5 WRITE(6,'(1X,A79)') LINE CLOSE(UNIT=1) ENDIF C*********************************************************************** C FIRST PRINT OUT HEADER FOR SIGHTING DISTANCE ESTIMATION AND C SOME EXPLANATORY INFORMATION C*********************************************************************** CALL HEADER (1) WRITE (6,160) CALL HEADER (2) C********************************************************************** C WRITE OUT ESTIMATORS CHOSEN C*********************************************************************** 10 WRITE (6,200) DO 20 I=1,NSD INDEX=SDEST(I) IF (INDEX.EQ.1) THEN WRITE (6,170) ELSE IF (INDEX.EQ.2) THEN WRITE (6,180) ELSE IF (INDEX.EQ.3) THEN WRITE (6,190) ENDIF 20 CONTINUE C*********************************************************************** C INITIALIZE AND SET UP SOME VARIABLES USED LATER ON C*********************************************************************** CZERO=0.0 VARC=0.0 PAR(1)=0.0 XN=FLOAT(N) SUM1=0.0 SUM2=0.0 DO 30 I=1,N TEMP=MESANG(I) IF (STATUS.EQ.1) TEMP=COMANG(I) SUM1=SUM1+TEMP SUM2=SUM2+TEMP*TEMP 30 CONTINUE THETM=SUM1/XN VART=(SUM2-(SUM1**2.)/XN)/((XN-1.)*XN) IF (WARN) VARN=XN C*********************************************************************** C CALCULATE THE MEAN OF THE RECIPROCAL SIGHTING DISTANCES AND ITS C VARIANCE. C*********************************************************************** RMEAN=0.0 VARR=0.0 DO 40 I=1,N 40 RMEAN=RMEAN+1./(SD(I)*XN) DO 50 I=1,N 50 VARR=VARR+(1./SD(I)-RMEAN)**2 VARR=VARR/(XN*(XN-1.)) C*********************************************************************** C CALCULATE ESTIMATES OF THE PARAMETERS AND DENSITY FOR EACH C ESTIMATOR REQUESTED C*********************************************************************** DO 150 I=1,NSD CALL HEADER (1) INDEX=SDEST(I) C*********************************************************************** C GENERALIZED HAYNE ESTIMATOR C*********************************************************************** IF (INDEX.EQ.1) THEN WRITE (6,270) WRITE (6,240) WRITE (6,280) CALL ELLIPS(SINA,CZERO,VARC,N) NPAR=1 C*********************************************************************** C MODIFIED HAYNE ESTIMATOR C*********************************************************************** ELSE IF (INDEX.EQ.2) THEN WRITE (6,270) WRITE (6,250) WRITE (6,280) IF ((THETM.GE.32.7).AND.(THETM.LE.45.0)) GO TO 70 WRITE (6,350) THETM GO TO 150 70 DELTA=(THETM-32.7)*.0813 CZERO=1.-.36338022*DELTA VARC=.0008726*VART NPAR=1 C*********************************************************************** C HAYNE ESTIMATOR C*********************************************************************** ELSE WRITE (6,270) WRITE (6,260) WRITE (6,280) NPAR=0 CZERO=1.0 VARC=0.0 ENDIF C*********************************************************************** C NEXT CALCULATE THE ESTIMATE OF F(0) AND STORE THE PARAMETER FOR TH C CALL TO PARAM WHICH WRITES OUT ALL OF THE PARAMETERS AND ASSOCIATE C STATISTICS C*********************************************************************** 90 PAR(1)=CZERO NPDI=NPD+I NPSET(NPDI)=NPAR VCMAT(1,1)=VARC FZERO=CZERO*RMEAN CVC=VARC/(CZERO*CZERO) CVR=VARR/(RMEAN*RMEAN) VARF=(FZERO*FZERO)*(CVR+CVC) CALL HEADER (2) CALL PARAM (XN,DNSTY,VARD,DLL,DUL) D(NPD+I)=DNSTY DL(NPD+I)=DLL DU(NPD+I)=DUL STD(NPD+I)=VARD C*********************************************************************** C NEXT OUTPUT THE PERCENTAGE OF THE VARIANCE OF F(0) WHICH IS C ATTRIBUTABLE TO ESTIMATION OF CZERO C*********************************************************************** IF (INDEX.EQ.3) GO TO 120 RATIO=100.*CVC/(CVC+CVR) WRITE (6,290) CVC,CVR,RATIO RATIO=100.*CVC/(CVC+CVR+(VARN/(XN*XN))) WRITE (6,390) RATIO KC=IFIX(SQRT(XN)) KC=MIN0(15,KC) DO 100 J=1,KC 100 SCUT(J)=FLOAT(J)/FLOAT(KC) IF (INDEX.EQ.2) GO TO 150 C*********************************************************************** C TEST ASSUMPTION OF THE GEN. HAYNE ESTIMATOR C*********************************************************************** CALL HEADER (1) WRITE (6,360) CALL SORT (U,N) WRITE (6,370) WRITE (6,450) (U(J),J=1,N) CALL HIST (U,N,SCUT,KC,FR) CHSQ=0.0 EXPECT=XN/FLOAT(KC) DO 110 J=1,KC SQ=FR(J)-EXPECT SQ=(SQ*SQ)/EXPECT 110 CHSQ=CHSQ+SQ WRITE (6,420) KC,CHSQ CALL CHI (KC-2,CHSQ,PROB) PROB=1.0-PROB KC2=KC-2 WRITE (6,430) PROB,KC2 C*********************************************************************** C PLOT THE CDF OF THE TRANSFORMED SINES AGAINST THE UNIFORM CDF C*********************************************************************** CALL HEADER (1) WRITE (6,380) CALL UTEST (U,N,N,2) WRITE (6,400) GO TO 150 C*********************************************************************** C TEST ASSUMPTIONS OF THE HAYNE ESTIMATOR THAT THE SIGHTING ANGLE C SINES ARE FROM A UNIFORM DISTRIBUTION. THREE TESTS AND A PLOT ARE C GIVEN. C*********************************************************************** 120 CALL HEADER (1) WRITE (6,210) WRITE (6,220) N IF (STATUS.NE.2) GO TO 130 WRITE (6,440) GO TO 150 C*********************************************************************** C TEST AVG. SINE = .5 C*********************************************************************** 130 CALL UTEST (SINA,N,N,1) C*********************************************************************** C TEST AVG. ANGLE = 32.7 C*********************************************************************** Z=(THETM-32.7)*SQRT(XN)/21.560218 WRITE (6,310) THETM,Z Z=Z*Z CALL CHI (1,Z,PROB) PROB=1.-PROB WRITE (6,320) PROB C*********************************************************************** C TEST IF THE SINES ARE FROM A UNIFORM DISTRIBUTION C BY THE CHI-SQUARE GOODNESS OF FIT TEST C*********************************************************************** CALL HIST (SINA,N,SCUT,KC,FR) EXPECT=XN/FLOAT(KC) CHSQ=0.0 DO 140 J=1,KC SQ=FR(J)-EXPECT SQ=(SQ*SQ)/EXPECT 140 CHSQ=CHSQ+SQ WRITE (6,330) KC,CHSQ CALL CHI (KC-1,CHSQ,PROB) PROB=1.-PROB KC1=KC-1 WRITE (6,340) PROB,KC1 C*********************************************************************** C PLOT CDF OF THE SINES AGAINST THE UNIFORM CDF C*********************************************************************** CALL HEADER (1) WRITE (6,300) CALL UTEST (SINA,N,N,2) WRITE (6,410) 150 CONTINUE RETURN C*********************************************************************** C FORMAT STATEMENTS C*********************************************************************** C 160 FORMAT (//10X,60('*')/10X,'*',58X,'*'/10X,'*',3X, 1 'Density Estimation from Sighting Distances and Angles',2X,'*'/ 2 ,10X,'*',58X,'*'/10X,60('*')) 170 FORMAT (/6X,'Generalized Hayne(GHYN)') 180 FORMAT (/6X,'Modified Hayne(MHYN)') 190 FORMAT (/6X,'Hayne(HAYN)') 200 FORMAT (//6X,'The estimators chosen for density estimation'/ 1 6X,'from sighting distances-angles are:') 210 FORMAT (//'0The Hayne estimator assumes that the flushing', 1' envelope is a circle. On the'/ 2' basis of this assumption, the distribution of the', 3' sighting angle is the cosine'/ 4' of the angle and the distribution of the', 5' sighting angle sine is uniform.'/ 6' Three tests are performed to test this assumption:', 7' 1) average sighting angle'/ 8' sine = 0.5 (the expected value under uniform', 9' distribution assumption), 2) aver-'/ a' age angle = 32.7 (expected value under cosine', B' assumption), and 3) chi-square'/ C' test that the PDF of sighting angle sines is uniform. Also,', D' a plot of the'/ E' sample CDF of sighting angle sines and the theoretical uniform', F' CDF is given.') 220 FORMAT (' The first two tests provide information about the', 1' central tendency of the'/ 2' sample distribution, whereas the third test examines the', 3' entire distribution.'/ 4' It is possible that the first two tests may not be significant', 5' and the chi-'/ 6' square test significant. The opposite is less likely to occur', 7' but is possible'/ 8' because the chi-square test is not very powerful.'// 9' Sample size = ',I4/'0Ordered Sighting Angle Sines') 240 FORMAT (21X,'*',5X,'Generalized Hayne Estimator',6X,'*') 250 FORMAT (21X,'*',7X,'Modified Hayne Estimator',7X,'*') 260 FORMAT (21X,'*',11X,'Hayne estimator ',11X,'*') 270 FORMAT (//21X,40('*')/21X,'*',38X,'*') 280 FORMAT (21X,'*',38X,'*'/21X,40('*')) 290 FORMAT ('0When the parameter A(1) (i.e., C) is estimated in a', 1' sighting distance model'/ 2' there is a corresponding increase in the variance. There', 3' is a trade-off'/ 4' between robustness and variance. The variance attributable', 5' to the estimation'/ 4' of the parameter is given below.'// 5 2X,'Squared Coefficient of Variation of A(1) (i.e., C) =',G11.4, 6/2X,'Squared Coefficient of Variation of the mean SD', 7' reciprocal (1/SD) =',G11.4// 7 2X,'Sampling variance of F(0) attributable to', 8' estimating A(1) is ',F5.2,'%.') 300 FORMAT (/9X,'Sample CDF of Sighting Angle Sines (*) and the', 1' Uniform CDF (u)') 310 FORMAT ('0Test of the hypothesis that the average sighting', 1' angle = 32.7 using the'/ 2' theoretical distribution of the sighting angle, which is', 3' the cosine of the'/ 4' angle. Thus the theoretical standard error of the average', 5' angle is '/ 6' equal to sqrt(464.843/n)'// 7' Average angle = ',F6.2,4X,'z-test value = ',F8.4) 320 FORMAT (' Significance level = ',f9.5//) 330 FORMAT ('0Chi-square goodness of fit test. The histogram of', 1' the sighting angle sines'/ 2' for the second cut point set in the data description section', 3' corresponds to the'/ 4' intervals used for this test. There were ',I2, 5' intervals used.'/'0Chi-square value = ',G11.4) 340 FORMAT (' Significance level = ',f9.5,3X, 1' Degrees of freedom = ',I2) 350 FORMAT (/'0Modified Hayne estimator is designed for situations', 1' when the average angle'/ 2' is between 32.7 and 45.0. The average angle is = ',f6.2,'. ', 3' Thus estimation will'/ 4' not be performed.') 360 FORMAT (/'0The Generalized Hayne estimator assumes that the', 1' flushing envelope is an ellipse'/ 2' and the ratio of the major and minor axes is', 3' constant and equal to the'/ 4' parameter C. The sines of the angle are not', 5' necessarily from a uniform dis-'/ 6' tribution but they can be transformed to a uniform', 7' distribution on the basis of'/ 8' the estimated parameter C if the model fits. A chi-square', 9' goodness of fit test'/ a' is performed to determine if the sample of transformed sines', B' comes from a uni-'/ C' form distribution. Also, the sample CDF of the transformed', D' sighting angle sines'/ E' and the uniform CDF are plotted together to allow a visual', F' inspection of the'/ G' fit.') 370 FORMAT (///' Ordered Transformed Sighting Angle Sines'/) 380 FORMAT (/2X,'Sample CDF of the Transformed Sighting Angle', 1' Sines (*) and the Uniform CDF (u)') 390 FORMAT (2X,'Sampling variance of D attributable to', 1' estimating A(1) is ',F5.2,'%.') 400 FORMAT (/34X,'Transformed Sighting Angle Sines') 410 FORMAT (/30X,'Sighting Angle Sines') 420 FORMAT (//' Chi-square goodness of fit test that the', 1' transformed sines'/ 2' come from a uniform distribution.'/3X,I3,' intervals of equal', 3' width were used.'/'0Chi-square value = ',G11.4) 430 FORMAT (' Significance level = ',F9.5/ 1' Degrees of freedom = ',I3) 440 FORMAT ('0*** WARNING - This assumption cannot be tested'/ 1' because sighting angles were not collected.') 450 FORMAT (5X,10F7.4) 490 FORMAT (///) END