C*********************************************************************** C THIS SUBROUTINE PROVIDES THE EXPECTED PROBABILITIES(CELL) WITHIN C AN INTERVAL AND THE PARTIALS OF THESE WITH RESPECT TO THE PARA- C METERS(PCELL). THESE VALUES ARE USED TO CALCULATE THE LOG- C LIKELIHOOD EQUATIONS AND THEIR PARTIALS. C C SUBROUTINES CALLED: NCNI,CHECK C*********************************************************************** SUBROUTINE CELLS (ITIME) C*********************************************************************** C DECLARATIONS C*********************************************************************** INCLUDE 'PARMTR.INC' DOUBLE PRECISION E(600), Y(600), E0, Z0, VALUE DOUBLE PRECISION EXPW, OEXPW, XX DOUBLE PRECISION PRF0A, PRF0B, RF0, CELL, PCELL INTEGER CNT LOGICAL GRP, POOL, PEST, SEST, DESC, DEF, CUTP, TRUNC, WARN LOGICAL HELP C*********************************************************************** C COMMON STATEMENTS C*********************************************************************** COMMON /NUM/ XL(MAXLIN), WIDTH, N, CNT, CONV(3), VARN, IDF, WARN COMMON /TSINTG/ E, Y COMMON /PART/ PCELL(MAXCEL,MAXPAR), CELL(MAXCEL) COMMON /INTER/ K, CUT(MAXCEL), F(MAXCEL), NCUT, NKC(5), 1 NK(5), RFREQ(5,MAXCEL), CCUT(5,MAXCEL) COMMON /OPTION/ GRP, POOL, PEST, SEST, DESC, DEF, CUTP, TRUNC, 1 HELP COMMON /INTEG/ RF0, PRF0A, PRF0B DOUBLE PRECISION PAR, VCMAT, G, XLL COMMON /DPAR/ PAR(MAXPR2), VCMAT(MAXPAR,MAXPAR), G(MAXPAR), 1 XLL, NPAR, INDEX COMMON /ERROR/ IER(4), LER LOGICAL LER C*********************************************************************** C DATA STATEMENTS C*********************************************************************** DATA PI /3.1415926536/ C*********************************************************************** C INITIALIZE VARIABLES C*********************************************************************** MNUM=120 RF0=0.0 PRF0A=0.0 PRF0B=0.0 CUTLEF=0.0 C*********************************************************************** C FOURIER SERIES ESTIMATOR C*********************************************************************** IF (INDEX.EQ.1) THEN DO 20 I=1,K XX=0.0 IF(NPAR.GT.0) THEN DO 10 J=1,NPAR XJ=FLOAT(J) A=SIN(XJ*PI*CUT(I)/WIDTH) B=SIN(XJ*PI*CUTLEF/WIDTH) PCELL(I,J)=((A-B)*WIDTH)/(XJ*PI) 10 XX=XX+PAR(J)*PCELL(I,J) ENDIF CELL(I)=((CUT(I)-CUTLEF)/WIDTH)+XX 20 CUTLEF=CUT(I) C*********************************************************************** C EXPONENTIAL POWER SERIES C*********************************************************************** ELSE IF (INDEX.EQ.2) THEN DO 70 I=1,K RWID=CUT(I)-CUTLEF DELTA=RWID/FLOAT(MNUM) DO 40 J=1,MNUM XX=CUTLEF+FLOAT(J)*DELTA XX=(XX/PAR(1))**PAR(2) Y(J)=XX 40 E(J)=DEXP(MAX(-460.D0,MIN(300.D0,-XX))) Z0=(CUTLEF/PAR(1))**PAR(2) E0=DEXP(MAX(-460.D0,MIN(300.D0,-Z0))) CALL NCNI (MNUM,E,E0,RWID,VALUE,IDIG) IF (IDIG.LT.6) IER(1)=-1 CELL(I)=VALUE RF0=RF0+CELL(I) IF (ITIME.EQ.1) THEN DO 50 J=1,MNUM 50 E(J)=Y(J)*E(J) E0=E0*Z0 CALL NCNI (MNUM,E,E0,RWID,VALUE,IDIG) IF (IDIG.LT.6) IER(1)=-1 PCELL(I,1)=VALUE PCELL(I,1)=PCELL(I,1)*PAR(2)/PAR(1) DO 60 J=1,MNUM 60 E(J)=E(J)*DLOG(MAX(1.D-200,Y(J))) IF (Z0.NE.0.0) Z0=DLOG(MAX(1.D-200,Z0)) E0=E0*Z0 CALL NCNI (MNUM,E,E0,RWID,VALUE,IDIG) IF (IDIG.LT.6) IER(1)=-1 PCELL(I,2)=-VALUE/PAR(2) PRF0A=PRF0A+PCELL(I,1) PRF0B=PRF0B+PCELL(I,2) ENDIF CUTLEF=CUT(I) 70 CONTINUE DO 80 I=1,K CELL(I)=CELL(I)/RF0 IF (ITIME.EQ.1) THEN PCELL(I,1)=(PCELL(I,1)-CELL(I)*PRF0A)/RF0 PCELL(I,2)=(PCELL(I,2)-CELL(I)*PRF0B)/RF0 ENDIF 80 CONTINUE C*********************************************************************** C EXPONENTIAL POLYNOMIAL ESTIMATOR C*********************************************************************** ELSE IF (INDEX.EQ.3) THEN K1=K OEXPW=1.0 IF (.NOT.TRUNC) K1=K-1 DO 110 I=1,K1 RWID=CUT(I)-CUTLEF DELTA=RWID/FLOAT(MNUM) E0=DEXP(MAX(-460.D0,MIN(300.D0, 1 (-PAR(1)-PAR(2)*CUTLEF)*CUTLEF))) DO 100 J=1,MNUM XX=CUTLEF+DELTA*FLOAT(J) 100 E(J)=DEXP(MAX(-460.D0,MIN(300.D0,(-PAR(1)-PAR(2)*XX)*XX))) CALL NCNI (MNUM,E,E0,RWID,VALUE,IDIG) IF (IDIG.LT.6) IER(1)=-1 CELL(I)=VALUE RF0=RF0+CELL(I) IF (ITIME.EQ.1) THEN EXPW=DEXP(MAX(-460.D0,MIN(300.D0, 1 (-PAR(1)-PAR(2)*CUT(I))*CUT(I)))) PCELL(I,1)=(EXPW-OEXPW+PAR(1)*CELL(I))/(2.D0*PAR(2)) PCELL(I,2)=(CUT(I)*EXPW-CUTLEF*OEXPW-CELL(I)- 1 PAR(1)*PCELL(I,1))/(2.D0*PAR(2)) PRF0A=PRF0A+PCELL(I,1) PRF0B=PRF0B+PCELL(I,2) OEXPW=EXPW ENDIF CUTLEF=CUT(I) 110 CONTINUE IF (.NOT. TRUNC) THEN CELL(K)=.5*SQRT(PI)*DEXP(MAX(-460.D0, 1 MIN(300.D0,PAR(1)*PAR(1)/(4.D0*PAR(2)))))/DSQRT(PAR(2)) RWID=PAR(1)/(2.D0*PAR(2)) DELTA=RWID/FLOAT(MNUM) E0=DEXP(MAX(-460.D0,MIN(300.D0,PAR(1)*PAR(1)/(4.D0*PAR(2))))) DO 120 J=1,MNUM XX=-PAR(1)/(2.D0*PAR(2))+FLOAT(J)*DELTA 120 E(J)=DEXP(MAX(-460.D0,MIN(300.D0,(-PAR(1)-PAR(2)*XX)*XX))) CALL NCNI (MNUM,E,E0,RWID,VALUE,IDIG) IF (IDIG.LT.6) IER(1)=-1 CELL(K)=CELL(K)-VALUE-RF0 RF0=RF0+CELL(K) IF (ITIME.EQ.1) THEN PCELL(K,1)=(-EXPW+PAR(1)*CELL(K))/(2.D0*PAR(2)) PCELL(K,2)=(-CUTLEF*EXPW-CELL(K)- 1 PAR(1)*PCELL(K,1))/(2.D0*PAR(2)) PRF0A=PRF0A+PCELL(K,1) PRF0B=PRF0B+PCELL(K,2) ENDIF ENDIF DO 140 I=1,K CELL(I)=CELL(I)/RF0 IF (ITIME.EQ.1) THEN PCELL(I,1)=(PCELL(I,1)-CELL(I)*PRF0A)/RF0 PCELL(I,2)=(PCELL(I,2)-CELL(I)*PRF0B)/RF0 ENDIF 140 CONTINUE C*********************************************************************** C NEGATIVE EXPONENTIAL ESTIMATOR C*********************************************************************** ELSE IF (INDEX.EQ.4) THEN K1=K OEXPW=1.0 IF (.NOT.TRUNC) K1=K-1 DO 160 I=1,K1 EXPW=DEXP(MAX(-460.D0,MIN(300.D0,-PAR(1)*CUT(I)))) CELL(I)=(OEXPW-EXPW)/PAR(1) RF0=RF0+CELL(I) IF (ITIME.EQ.1) THEN PCELL(I,1)=(CUT(I)*EXPW-CUTLEF*OEXPW-CELL(I))/PAR(1) PRF0A=PRF0A+PCELL(I,1) ENDIF OEXPW=EXPW CUTLEF=CUT(I) 160 CONTINUE IF (.NOT. TRUNC) THEN CELL(K)=OEXPW/PAR(1) RF0=RF0+CELL(K) IF (ITIME.EQ.1) THEN PCELL(K,1)=(-CUTLEF*OEXPW-CELL(K))/PAR(1) PRF0A=PRF0A+PCELL(K,1) ENDIF ENDIF DO 180 I=1,K CELL(I)=CELL(I)/RF0 IF (ITIME.EQ.1) THEN PCELL(I,1)=(PCELL(I,1)-CELL(I)*PRF0A)/RF0 ENDIF 180 CONTINUE C*********************************************************************** C HALF-NORMAL ESTIMATOR C*********************************************************************** ELSE IF (INDEX.EQ.5) THEN K1=K OEXPW=1.0 IF (.NOT.TRUNC) K1=K-1 DO 210 I=1,K1 RWID=CUT(I)-CUTLEF DELTA=RWID/FLOAT(MNUM) E0=DEXP(MAX(-460.D0,MIN(300.D0,-PAR(1)*CUTLEF*CUTLEF))) DO 200 J=1,MNUM XX=CUTLEF+DELTA*FLOAT(J) 200 E(J)=DEXP(MAX(-460.D0,MIN(300.D0,-PAR(1)*XX*XX))) CALL NCNI (MNUM,E,E0,RWID,VALUE,IDIG) IF (IDIG.LT.6) IER(1)=-1 CELL(I)=VALUE RF0=RF0+CELL(I) IF (ITIME.EQ.1) THEN EXPW=DEXP(MAX(-460.D0,MIN(300.D0,-PAR(1)*CUT(I)*CUT(I)))) PCELL(I,1)=(CUT(I)*EXPW-CUTLEF*OEXPW-CELL(I))/(2.D0*PAR(1)) PRF0A=PRF0A+PCELL(I,1) OEXPW=EXPW ENDIF CUTLEF=CUT(I) 210 CONTINUE IF (TRUNC) GO TO 220 CELL(K)=SQRT(PI)*.5/DSQRT(PAR(1))-RF0 RF0=RF0+CELL(K) IF (ITIME.EQ.1) THEN PCELL(K,1)=(-CUTLEF*OEXPW-CELL(K))/(2.D0*PAR(1)) PRF0A=PRF0A+PCELL(K,1) ENDIF 220 DO 230 I=1,K CELL(I)=CELL(I)/RF0 IF (ITIME.EQ.1) THEN PCELL(I,1)=(PCELL(I,1)-CELL(I)*PRF0A)/RF0 ENDIF 230 CONTINUE ENDIF CALL CHECK (CELL,K) XLL=0.0 DO 240 I=1,K 240 XLL=XLL+F(I)*DLOG(MAX(1.D-200,CELL(I))) IF (ITIME.EQ.1) THEN C C FIRST COMPUTE G(I),I=1,...NPAR,WHICH ARE THE LOG-LIKELIHOOD EQUATION C DO 250 J=1,NPAR G(J)=0.0 DO 260 I=1,K 260 G(J)=G(J)+(F(I)*PCELL(I,J)/CELL(I)) 250 CONTINUE C C NOW COMPUTE THE INFORMATION MATRIX C DO 270 J=1,NPAR DO 270 M=J,NPAR VCMAT(J,M)=0.0 DO 280 I=1,K 280 VCMAT(J,M)=VCMAT(J,M)+PCELL(I,J)*PCELL(I,M)/CELL(I) VCMAT(J,M)=DBLE(CNT)*VCMAT(J,M) VCMAT(M,J)=VCMAT(J,M) 270 CONTINUE ENDIF RETURN END