SUBROUTINE MODLB2(TITLE,BCC) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER IDIM PARAMETER (IDIM=50) INTEGER NITER,I,J,K,L,S,MAXITR,YR(IDIM),BCC,IPRNT,ITRY CHARACTER*128 TITLE DIMENSION Q(IDIM),THETA(IDIM+2),THETA1(IDIM+2),H(IDIM+2,IDIM+2), , G(IDIM+2),VCP(IDIM,IDIM),COVA(IDIM) COMMON/BLK1/M(IDIM),N(IDIM),NY(IDIM),B(IDIM),R(IDIM),RY(IDIM), , NS(IDIM),SY(IDIM),YR COMMON /BLK3/S,FIMA(2),FYMA(2),EPSLON,MAXITR,IPRNT COMMON /BLK4/XA(IDIM),XY(IDIM),P(IDIM),XAD(IDIM),XYD(IDIM),PD(1) C C THIS SUBROUTINE COMPUTES ESTIMATES OF ADULT & YOUNG C SURVIVAL RATE (PHI & PHI') ASSUMING THEY ARE CONSTANT C AND TIME SPECIFIC CAPTURE Probability. C BCC=0 WRITE(7,10)CHAR(12),TITLE 10 FORMAT(A1,A128/ //' Model B2 - Constant survival rate per unit time, ', 2'time-specific capture probabilities.'/ 3/' Parameters: PHI = Adult survival rate per unit time'/ 5 12X,6H PHI',' = Young survival rate per unit time'/ 6 ' P(i) = Adult capture probability in sample i'/) ITRY=0 7 ITRY=ITRY+1 C C FIRST, COMPUTE STARTING VALUE OF THETA C XA(S)=.1D1 AVGP=.0D0 DO 20 I=S-2,1,-1 J=I+1 P(J)=.5D0 IF(R(J)*(M(J)+B(J)).GT..0D0)P(J)=M(J)/(M(J)+NS(J)*B(J)/R(J)) IF(P(J).GE..9)P(J)=.9 20 AVGP=AVGP+P(J) P(S)=AVGP/DBLE(FLOAT(S-2)) THETA(1)=FIMA(ITRY) DO 30 I=S-1,1,-1 J=I+1 Q(J)=.1D1-P(J) XA(I)=.1D1-(FIMA(ITRY)*(.1D1-Q(J)*XA(J))) XY(I)=.1D1-(FYMA(ITRY)*(.1D1-Q(J)*XA(J))) 30 THETA(J)=XA(I) THETA(S+1)=FYMA(ITRY) FI=FIMA(ITRY) FY=FYMA(ITRY) C C WRITE STARTING VALUES C WRITE(7,40)THETA(1),THETA(S+1) 40 FORMAT(/' Starting values of PHI (Adult, Young):',2F10.4) WRITE(7,50)(P(I),I=2,S) 50 FORMAT(/' Starting values of p:'/(10(1X,F8.4))) NITER=0 GO TO 80 C C START OF ITERATIVE PROCEDURE C SET SURV & CAPT. PROBS = THETA C 60 FI=THETA(1) FY=THETA(S+1) DO 70 I=S,2,-1 XA(I-1)=THETA(I) P(I)=((.1D1-XA(I-1))/FI-(.1D1-XA(I)))/XA(I) Q(I)=.1D1-P(I) XY(I-1)=.1D1-FY*(.1D1-Q(I)*XA(I)) 70 CONTINUE 80 CONTINUE C C COMPUTE GRADIANT VECTOR (G) & INVERSE OF VAR-COV MATRIX(H) C SUM=.0D0 SUM2=.0D0 SUM3=RY(1)-(.1D1-XY(1))/XY(1)*(SY(1)-RY(1)) SUM4=(.1D1-XY(1))*SY(1)/XY(1) DO 90 I=2,S DO 90 J=2,S 90 H(J,I)=.0D0 DO 100 I=2,S-1 SUM=SUM+(B(I)/Q(I)-M(I)*(.1D1-XA(I))/P(I))/XA(I) SUM2=SUM2+(SY(I)*(.1D1-XY(I))/XY(I)+(.1D1-Q(I)*XA(I))/ / XA(I)/XA(I)*(M(I)*(.1D1-XA(I))/P(I)/P(I)+ + B(I)/Q(I)/Q(I))) SUM3=SUM3+RY(I)-(.1D1-XY(I))/XY(I)*(SY(I)-RY(I)) SUM4=SUM4+(.1D1-XY(I))*SY(I)/XY(I) 100 CONTINUE C P(1)=.5D0 DO 110 I=1,S-1 J=I+1 BQ=.0D0 BQQ=.0D0 IF(B(J).NE..0D0)BQ=B(J)/Q(J) IF(B(J).NE..0D0)BQQ=BQ/Q(J) G(J)=(SY(I)-RY(I))/XY(I)*FY/FI+(M(I)*Q(I)/P(I)-B(I)+ + NS(I)-R(I))/XA(I)-(M(J)/P(J)-BQ)/FI/XA(J) IF(I.GT.1)H(I,J)=-M(I)/FI/P(I)/P(I)/XA(I)/XA(I) H(J,I)=H(I,J) H(1,J)=(FY/FI*SY(I)/XY(I)-(.1D1-Q(I)*XA(I))*M(I)/P(I)/P(I)/ / XA(I)/XA(I)+(M(J)*(.1D1-XA(J))/P(J)/P(J)+ + BQQ)/FI/XA(J)/XA(J))/FI H(J,1)=H(1,J) H(J,J)=FY*FY/FI/FI*SY(I)/XY(I)+NS(I)/XA(I)+ + (M(I)/P(I)/P(I)-M(I)-B(I))/XA(I)/XA(I)+ + (M(J)/P(J)/P(J)+BQQ)/FI/FI/XA(J)/XA(J) IF(I.GT.1)H(I,S+1)=-SY(I-1)/FI/XY(I-1) 110 H(S+1,I)=H(I,S+1) G(1)=SUM/FI G(S+1)=SUM3/FY H(1,1)=SUM2/FI/FI+SY(1)*(.1D1-XY(1))/FI/FI/XY(1) H(1,S+1)=-SUM4/FI/FY H(S+1,1)=H(1,S+1) H(S+1,S+1)=SUM4/FY/FY H(S,S+1)=-SY(S-1)/FI/XY(S-1) H(S+1,S)=H(S,S+1) C C COMPUTE NEW THETA = THETA + INV(H)*G C CALL MATINV(S+1,H,IDIM+2) CALL MMULT(H,G,THETA1,S+1,S+1,1,IDIM+2,IDIM+2,1) CALL MATADD(THETA1,THETA,THETA1,S+1) C C CHECK FOR CONVERGENCE C DO 120 I=1,S+1 IF(DABS(THETA(I)-THETA1(I)).GT.EPSLON)GO TO 130 120 CONTINUE GO TO 180 C C NOT CONVERGED... INCREASE # ITERATIONS & SET THETA=NEW THETA C 130 NITER=NITER+1 IF(NITER.GT.MAXITR)GO TO 150 IF(THETA1(1).LT..0D0.OR.THETA1(S+1).LT..0D0)GO TO 150 DO 140 I=1,S+1 IF(DABS(THETA1(I)).GT..1D2)GO TO 150 140 THETA(I)=THETA1(I) GO TO 60 C C CONVERGENCE FAILURE... OUTPUT MESSAGE & QUIT C 150 WRITE(7,160)NITER,(THETA(I),I=1,S+1) 160 FORMAT(/' Convergence failure after ',I6,' iterations'/ /' THETA = '/(10F12.6)) WRITE(7,170)(THETA1(I),I=1,S+1) 170 FORMAT(/' THETA1 = '/(10F12.6)) IF(ITRY.LE.1.AND.FIMA(2).GT..0D0)GO TO 7 RETURN C C CONVERGED... OUTPUT RESULTS C 180 BCC=1 SE=XSQRT(H(1,1)) C1=FI-SE*.196D1 C2=FI+SE*.196D1 WRITE(7,190)NITER,FI,H(1,1),SE,C1,C2,H(1,S+1) 190 FORMAT(/' Final values after ',I6,' iterations'// /' Parameter Estimate Variance Std.Error', ,' 95% confidence interval Covariance w/ PHI',3X, ,'Covariance W/ PHI',1H'/1X,110('-')/ /5X,'PHI',4X,F8.4,F10.4,F11.4,F10.4,' - ',F8.4,31X,F13.8/) SE=XSQRT(H(S+1,S+1)) C1=FY-SE*.196D1 C2=FY+SE*.196D1 WRITE(7,200)FY,H(S+1,S+1),SE,C1,C2,H(1,S+1) 200 FORMAT(5X,4HPHI',3X,F8.4,F10.4,F11.4,F10.4,' - ',F8.4,10X,F13.8/) AVGP=.0D0 AVGVAR=.0D0 DO 240 I=2,S DO 210 J=I,S T1=(.1D1-Q(I)*XA(I))/FI T2=(.1D1-Q(J)*XA(J))/FI VCP(I,J)=(T1*(T2*H(1,1)+H(1,J)/FI)+ + (T2*H(1,I)+H(I,J)/FI)/FI)/XA(I)/XA(J) IF(I.LT.S)VCP(I,J)=VCP(I,J)-Q(I)/XA(I)/XA(J)* * (T2*H(1,I+1)+H(I+1,J)/FI) IF(J.LT.S)VCP(I,J)=VCP(I,J)-Q(J)/XA(I)/XA(J)* * (T1*H(1,J+1)+H(I,J+1)/FI) IF(J.LT.S.AND.I.LT.S)VCP(I,J)=VCP(I,J)+Q(I)*Q(J)/XA(I)/XA(J)* * H(I+1,J+1) AVGVAR=AVGVAR+VCP(I,J) 210 CONTINUE AVGP=AVGP+P(I) SE=XSQRT(VCP(I,I)) C1=P(I)-SE*.196D1 C2=P(I)+SE*.196D1 COVA(I)=-(H(1,I)+(.1D1-Q(I)*XA(I))*H(1,1))/FI/XA(I) IF(I.LT.S)COVA(I)=COVA(I)+Q(I)*H(1,I+1)/XA(I) COVY=-(H(I,S+1)+(.1D1-Q(I)*XA(I))*H(1,S+1))/FI/XA(I) IF(I.LT.S)COVY=COVY+Q(I)*H(I+1,S+1)/XA(I) WRITE(7,220)YR(I),P(I),VCP(I,I),SE,C1,C2,COVA(I),COVY 220 FORMAT(5X,'p(',I4,')',F8.4,F10.4,F11.4,F10.4,' - ', , F8.4,10X,F13.8,8X,F13.8) IF(MOD(YR(I),5).EQ.0)WRITE(7,230) 230 FORMAT(I5) 240 CONTINUE AVGP=AVGP/(S-1) AVGVAR=AVGVAR/(S-1) SE=XSQRT(AVGVAR) C1=AVGP-SE*.196D1 C2=AVGP+SE*.196D1 WRITE(7,250)AVGP,AVGVAR,SE,C1,C2 250 FORMAT(1X,62('-')/' Avg p',6X,F8.4,F10.4,F11.4,F10.4,' - ', , F8.4/1X,62('-')) CALL B2MNB(TITLE,FI,FY,P,XA,XY,H,VCP,COVA) RETURN END