C$CONTROL USLINIT,SEGMENT 'MAIN' C PROGRAM JOLLY C Changes: 01/22/91 Fixed output of Covariance matrix when nyrs>10 in mod A C 08/28/90 Restricted between-period resighting data to model A C 08/01/90 Increased dimensions to 50. C 03/20/90 Fixed bug in pooling routine. C 01/31/90 Fixed bug in GOFM2 when N2PTR=0 C 02/24/89 Added FREE format to input routine C 10/07/88 Added codes 3 & 4 to capture-history record C format. C 7/15/88 Added printout of p,SE(p) to model A' C 7/07/88 Fixed bug in model 2 covariance of PHI**1/T(i) C 5/10/88 Modified A' gof test : 3 -> 4 cols. C Added test between models A & A' C 4/28/88 Fixed MOD A' var(n) to not include terms < 0. C 1/19/88 Fixed chi-square routine for high chi-sq values C 12/01/87 Changed input routine to get title,nyrs... C from input data. C 7/10/87 Added menu driven input routine C 7/01/87 Added billions and billions of IF statements C to avoid divide by zero in models A, A', and 2. C 1/29/87 Added PHI**1/T(I) & var. to mods A, A' & 2 C 12/16/86 Converted to HP FORTRAN 77. C Fixed model A output of phi**t(i) C Inserted routine in model B to try diff strt vals C 8/11/86 Variable T(IDIM) changed to double precision, C PHI(I)**1/T(I) printed out in model A. C 4/2/86 B-Table GOF test for model A' modified to compute C obs & expected values of # unmarked C 2/18/86 Old gof test fixed where R(J) or M(J)=0, C Old gof test used if 1st gof test DF=0, C IDFA=0 added to beg. of subr GOF2. C 1/22/86 var MODA passed to subr MODELA to enable C 2nd GOF test if MODA=2. C PROGRAM 'JOLLY' C THIS PROGRAM IS USED TO COMPUTE SURVIVAL RATE, POPULATION C SIZE, AND IMMIGRATION RATE ESTIMATES FROM CARTURE-RECAPTURE C DATA. THE PROGRAM COMPUTES THESE STATISTICS USING THE MODEL C SUGGESTED BY JOLLY (1965; BIOMETRIKA 52:225-246) WHICH C ALLOWS BOTH DEATH AND IMMIGRATION. A Goodness-of-fit TEST IS C ALSO PROVIDED TO DETERMINE IF THIS MODEL IS APPROPRIATE. A C SECOND MODEL IS AVAILABLE WHICH ALLOWS FOR DEATH BUT NO IMMIGRATION. C THE SAME STATISTICS (EXCEPT FOR THE ESTIMATE OF C IMMIGRATION), INCLUDING A Goodness-of-fit TEST, ARE COMPUTED IN C THIS MODEL. C TWO REDUCED PARAMETER MODELS (B & D) ARE ALSO AVAILABLE. MODEL C B ASSUMES CONSTANT SURVIVAL RATE AND TIME-SPECIFIC CAPTURE C PROBABILITY, AND MODEL D ASSUMES CONSTANT SURVIVAL RATE AND C CONSTANT CAPTURE PROBABILITY. C C IF THERE ARE ANY QUESTIONS OR PROBLEMS WHITH THIS PROGRAM CONTACT: C JAMES E. HINES OR JAMES D. NICHOLS C U. S. FISH & WILDLIFE SERVICE C PATUXENT WILDLIFE RESEARCH CENTER C LAUREL, MARYLAND 20708 C IMPLICIT DOUBLE PRECISION (A-Z) INTEGER I,J,K,L,M,NTIME,IDF,MAXITR,IDFA,IDFAP,PRNT1,PRNT2, ,IDF2,MODA,MOD1A,MOD2,MODB,MODD,YR,CHOICE,A1CC,M2CC,IDFA2,IDIM INTEGER*2 ARGC PARAMETER (IDIM=50) CHARACTER*3 YESNO(2),ANS*1,MDLS*20,ARGV*72 CHARACTER*40 FNAME,FNAME2,TMPCH*1 DIMENSION Z(IDIM),NNIJ(IDIM),STRT(IDIM),ZPRIME(IDIM),T(IDIM), , TITLE(16) COMMON /BLK0/NTIME,YR(IDIM),TITLE COMMON /BLK3/TABLE(IDIM,IDIM),NM(IDIM),NU(IDIM),NN(IDIM),NS(IDIM) DATA T,NNIJ,Z,ZPRIME,STRT/IDIM*.1D1,IDIM*.0D0,IDIM*.0D0, , IDIM*.0D0,IDIM*.0D0/ DATA IDFA,IDFAP,IDF2/3*0/ DATA CHIAP,CHI2,PRNT1,PRNT2/2*.0D0,2*2/ DATA MODA,MOD1A,MOD2,MODB,MODD/5*1/ DATA MDLS,FNAME,FNAME2/'A,B,D ','JOLLYINP ','JOLLYOUT'/ DATA YESNO/'YES','NO '/ IF(ARGC().GE.2)FNAME=ARGV(1) IF(ARGC().GE.3)FNAME2=ARGV(2) IF(ARGC().GE.2)GO TO 100 1 CONTINUE PRINT 5 5 FORMAT(/////17X,' Program JOLLY <01/24/91>'/ / 17X,' ',24('-')/) PRINT 6,FNAME,FNAME2,MDLS,YESNO(PRNT1),YESNO(PRNT2) 6 FORMAT(17X,' 1) Data filename:',A40/ / 17X,' 2) Name of print file:',A40/ / 17X,' 3) Models:',A20/17X,' 4) Print definitions:',A3/ /17X,' 5) Print Var-Cov of N:',A3/17X,' 6) Run...'/17X,' 7) Quit'/ //17X,' Choice(1...7)?') READ(*,*,END=101,ERR=900)CHOICE GO TO (60,90,40,50,55,100,800),CHOICE 8 PRINT *,'(',CHOICE,') Please enter a number between 1 and 7!' PRINT *,' ** Press to continue...' READ(*,FMT='(A1)',ERR=900)TMPCH GO TO 1 40 PRINT *,'Enter desired models (eg A,A'',2,B,D):' READ(*,FMT='(A20)',ERR=900)MDLS GO TO 1 50 PRNT1=3-PRNT1 GO TO 1 55 PRNT2=3-PRNT2 GO TO 1 60 PRINT *,'Enter new filename:' READ(*,FMT='(A40)',ERR=900)FNAME GO TO 1 C 90 PRINT *,'Name of print file?' READ(*,FMT='(A40)',ERR=900)FNAME2 GO TO 1 C CALL SELECTED MODEL SUBROUTINES 100 CONTINUE 101 MAXITR=100 EPSLON=.1D-4 OPEN(1,FILE=FNAME,STATUS='OLD',ERR=950) OPEN(7,FILE=FNAME2) IF(PRNT1.EQ.1)CALL DEFINE C IF(INDEX(MDLS,'ax')+INDEX(MDLS,'AX').GT.0)MDLS='AX' CALL BTABLE(TITLE,MDLS,NTIME,YR,T,NNIJ,Z,ZPRIME,STRT,PRNT1,PRNT2) C C COMPUTE ESTIMATES UNDER MODEL A (MODEL 1) C CALL MODELA(T,NM,NN,NS,Z,NNIJ,PHI,PEST,CHIA,IDFA,MODA,PRNT2, , CHIA2,IDFA2) IF(INDEX(MDLS,'a''')+INDEX(MDLS,'A''').GT.0)THEN CALL M1A(T,NU,NN,NS,NNIJ,ZPRIME,MOD1A,CHIAP,IDFAP,PRNT2) A1CC=1 ELSE A1CC=0 ENDIF IF(INDEX(MDLS,'2').GT.0)CALL MODEL2(T,NM,NN,NS,Z,CHI2,IDF2,M2CC) C IF(INDEX(MDLS,'d')+INDEX(MDLS,'D').GT.0) CALL MODELD X (T,NM,NN,NS,NNIJ,Z,STRT,EPSLON,MAXITR,PHI,PEST,MODD,PRNT2) C IF(INDEX(MDLS,'B')+INDEX(MDLS,'b').GT.0)CALL MODELB X (T,NM,NN,NS,NNIJ,Z,STRT,EPSLON,MAXITR,PHI,MODB,PRNT2) CALL TESTS(NM,NN,NS,NNIJ,Z,A1CC,M2CC,MODD,MODB,CHIA,IDFA, , CHIAP,IDFAP,CHI2,IDF2,CHIA2,IDFA2) 800 STOP 900 PRINT *,'Input error - press RETURN to continue:' 901 READ(*,FMT='(A1)')TMPCH GO TO 1 950 PRINT *,'Error: Input file not found! Press RETURN to continue:' IF(ARGC().LT.2)GO TO 901 STOP END C$CONTROL SEGMENT 'BTBLE' FUNCTION UPPER(STR) CHARACTER*250 STR,UPPER DO 10 I=1,LEN(STR) UPPER(I:I)=STR(I:I) IF(STR(I:I).GE.'a'.AND.STR(I:I).LE.'z') ) UPPER(I:I)=CHAR(ICHAR(STR(I:I))-ICHAR('a')+ICHAR('A')) 10 CONTINUE RETURN END C$CONTROL SEGMENT 'BTBLE' SUBROUTINE BTABLE(TITL,MDLS,NT,YR,T,R,ZZ,ZP,STRT,PRNT,PRNT2) C SUBROUTINE READS INDIVIDUAL CAPTURE HISTORY RECORDS AND C CREATES A B-TABLE FROM THIS DATA. A SECOND IS ALSO C CREATED WHICH INCLUDES ONLY RE-RECAPTURES. C------------------------------------------------------------------- C THE NUMBER OF FIRST CAPTURES WHICH WERE C RELEASED (V(I)) IS STORED ON THE DIAGONAL OF BTBLE ARRAY. C IMPLICIT INTEGER (A-Z) PARAMETER (IDIM=50) INTEGER*2 YEAR,MNTH,DAY,HR,MINUTE,SEC,HNDSEC DOUBLE PRECISION T(IDIM),DTABLE,CA,CM,CN,Z,CMC,NM,NU,NN,NS, ,TITL(16),ZC,R(IDIM),ZZ(IDIM),ZP(IDIM),STRT(IDIM),EPS CHARACTER*20 MDLS,FMTYPE(4),DASH(1)*5 CHARACTER*4 KEYWD(12),ACODE*1,YCODE*1,CAGE*1 CHARACTER*250 FMT,LINE,UPPER,TITLE*128 C C SEE SUBROUTINE GOF2 FOR DEFINITIONS OF THE FOLLOWING C VARIABLES. C INTEGER CAP(IDIM*2),YR(IDIM),NT COMMON /BLK3/DTABLE(IDIM,IDIM),NM(IDIM),NU(IDIM),NN(IDIM),NS(IDIM) COMMON /BLK2/CA(IDIM,2),CM(IDIM,2),CN(IDIM,2),Z(IDIM,2),CMC(IDIM), , ZC(IDIM) C *** INITIALIZE DATA IERR,IREC/0,0/ ,DASH/'-----'/ DATA FMTYPE/'CAPTURE HISTORY (1)','CAP/REL HISTORY (11)', , 'B-TABLE','TOTALS'/ C READ FORMAT OF INPUT DATA KEYWD/'FORM','FIRS','TYPE','ITER','CONV','INTE', , 'PERI','MODE','COVA','STAR','DEFI','FILE'/ C CALL GETDAT(YEAR,MNTH,DAY) CALL GETTIM(HR,MINUTE,SEC,HNDSEC) WRITE(TITLE,901)MNTH,DAY,YEAR-1900,HR,MINUTE 901 FORMAT(105X,2(I2,'/'),I2,2X,I2,':',I2) TITLE(1:24)=' JOLLY (Ver:01/24/91) ' YR(1)=1 NT=0 MAXITR=100 EPS=.0001 IFMT=1 FMT='****' FLAGAX=0 DO 101 I=1,IDIM NM(I)=.0D0 NU(I)=.0D0 NN(I)=.0D0 NS(I)=.0D0 DO 101 J=1,IDIM 101 DTABLE(J,I)=.0D0 DO 102 I=1,IDIM CMC(I)=.0D0 ZC(I)=.0D0 T(I)=.1D1 DO 102 J=1,2 CA(I,J)=.0D0 CM(I,J)=.0D0 CN(I,J)=.0D0 102 Z(I,J)=.0D0 42 READ(1,'(A)')LINE IREC=IREC+1 IF(PRNT.EQ.1)WRITE(7,43)LINE 43 FORMAT(' INPUT==>',A) IF(INDEX(UPPER(LINE),'TITL').GT.0)GO TO 2 KW=0 DO 45 I=1,12 IF(INDEX(UPPER(LINE),KEYWD(I)).GT.0.AND.KW.EQ.0)KW=I 45 CONTINUE LEQ=INDEX(LINE,'=')+1 GO TO (1,3,4,5,6,7,8,9,10,11,8,12),KW+1 1 BACKSPACE 1 GO TO 79 2 TITLE(24:103)=LINE(LEQ:) IF(PRNT.EQ.1)WRITE(7,*)'CONTROL RECORD:TITLE=',TITLE(24:103) GO TO 42 3 FMT=LINE(LEQ:) GO TO 42 4 YR(1)=INUM(LINE(LEQ:)) IF(PRNT.EQ.1)WRITE(7,*)'CONTROL RECORD:FIRST=',YR(1) GO TO 42 5 IF(INDEX(UPPER(LINE),'TABLE').GT.0)IFMT=3 IF(INDEX(UPPER(LINE),'TOTAL').GT.0)IFMT=4 IF(INDEX(UPPER(LINE),'REL').GT.0)IFMT=2 IF(PRNT.EQ.1)WRITE(7,*)'CONTROL RECORD:TYPE =',FMTYPE(IFMT) GO TO 42 6 MAXITR=INUM(LINE(LEQ:)) IF(PRNT.EQ.1)WRITE(7,*)'CONTROL RECORD:ITER =',MAXITR GO TO 42 7 EPS=.1D2**(-INUM(LINE(LEQ:))) IF(PRNT.EQ.1)WRITE(7,*)'CONTROL RECORD:CONV =',EPS GO TO 42 8 I1=INDEX(LINE,'(') I2=INDEX(LINE,')') I=1 IF(I2-I1-1.GT.0)I=INUM(LINE(I1+1:I2-1)) IF(KW.EQ.6)THEN CALL GETARY(LINE(LEQ:),T(I)) IF(PRNT.EQ.1)WRITE(7,*)'INPUT==>',(T(I),I=1,NT-1) ELSE CALL GETARY(LINE(LEQ:),STRT(I)) IF(PRNT.EQ.1)WRITE(7,*)'INPUT==>',(STRT(I),I=1,NT-1) ENDIF GO TO 42 9 NT=INUM(LINE(LEQ:)) IF(PRNT.EQ.1)WRITE(7,*)'CONTROL RECORD:PERIODS=',NT GO TO 42 10 MDLS=LINE(LEQ:) GO TO 42 11 PRNT2=1 GO TO 42 12 PRNT=1 GO TO 42 13 CLOSE(1) OPEN(1,FILE=LINE(LEQ:),STATUS='OLD') C 79 WRITE(7,78)CHAR(12),TITLE 78 FORMAT(A1,A128) WRITE(7,39)NT,FMT,MAXITR,EPS 39 FORMAT(/' Number of sampling periods = ',I5/ / ' Input format = ',A70/ /' Max iterations = ',I6/' Convergence criterion = ',E10.4) WRITE(7,204)(T(I),I=1,NT-1) 204 FORMAT(/' Time period interval lengths = '/(10(1X,F12.8))) WRITE(7,*)'Data type:',FMTYPE(IFMT),' - ',FMT IF(NT.GT.IDIM)THEN WRITE(7,*)'** Error: Number of sampling periods > ',IDIM WRITE(7,*)'** Recompile program with increased IDIM' STOP ENDIF C C 3 POSSIBLE WAYS TO ENTER DATA C 1) DEFAULT(CAPTURE HISTORY RECORDS) EACH RECORD CONTAINS STATUS C (CAPTURED,RELEASED) OF ANIMAL FOR EACH WEEK OF EXPERIMENT. C EG. 11 00 10 : CAP & REL IN TIME 1, CAP & NOT REL IN 3. C 2) CAPTURE-HISTORY RECORDS WITH 1 DIGIT PER TIME-PERIOD C EG. 1 0 2 : CAP & REL IN TIME 1, CAP & NOT REL IN 3. C 3) B-TABLE FORMAT- B-TABLE IS ENTERED DIRECTLY. C C **************************************************** C **** NEW CODE ADDED TO CAPTURE HISTORY RECORDS **** C **** CAP(I)=3 -> CAPTURED AFTER PERIOD I, BUT **** C **** BEFORE PERIOD I+1 **** C **************************************************** C 40 IF(IFMT.EQ.3)GO TO 140 IF(IFMT.EQ.4)GO TO 195 C *** DEFAULT FORMAT IF(FMT(1:1).NE.'(')THEN READ(1,*,END=200,ERR=900)(CAP(I),I=1,NT*IFMT),CNT ELSE READ(1,FMT,END=200,ERR=900)(CAP(I),I=1,NT*IFMT),CNT ENDIF IF(IFMT.EQ.2)THEN DO 44 I=1,2*NT,2 CAP(I/2+1)=CAP(I) 44 IF(CAP(I).EQ.1.AND.CAP(I+1).EQ.0)CAP(I/2+1)=2 ENDIF IREC=IREC+1 IF(CNT.LE.0)CNT=1 LCAP=-1 LCAP2=-1 FRST=999 DO 80 I=1,NT IF(CAP(I).GE.3)FLAGAX=1 IF(CAP(I).GT.0)LCAP2=I IF(MOD(CAP(I),3).GT.0)THEN LCAP=I IF(I.LT.FRST)FRST=I NN(I)=NN(I)+CNT IF(CAP(I).NE.2)NS(I)=NS(I)+CNT IP1=I+1 IF(IP1.LE.NT)THEN J=0 DO 70 JJ=IP1,NT 70 IF(MOD(CAP(JJ),3).GT.0.AND.J.EQ.0)J=JJ IF(J.GT.0)THEN DTABLE(J,I)=DTABLE(J,I)+CNT IF(I.NE.FRST)DTABLE(I,J)=DTABLE(I,J)+CNT ENDIF ENDIF ENDIF 80 CONTINUE C C *** COMPUTE STATS NEEDED TO COMPUTE CONTINGENCY TABLES C CA(I,1)=# MARKED WHICH WERE CAUGHT IN SAMPLE I & after C CA(I,2)=# UNMARKED " " " " " " " " C CN(I,1)=# MARKED IN I, & NOT CAUGHT LATER C CN(I,2)=# UNMARKED " " " " " " C CM(I,1)=# CAUGHT AS MARKED IN SAMPLE I & CAUGHT before C CM(I,2)=# " " " " " " " NOT " " C CMC(I) =# CAUGHT IN SAMPLES I & I-1 AND CAUGHT before I-1 C ZC(I) =# CAUGHT before I-1,IN I-1, NOT IN I, & after I C IF(LCAP.LT.FRST)GO TO 40 IF(INDEX(MDLS,'AX').GT.0.AND.CAP(LCAP2).GE.3)THEN LCAP2=LCAP2+1 ELSE LCAP2=LCAP ENDIF DO 130 I=FRST,NT IF(MOD(CAP(I),3).GT.0)THEN J=1 IF(I.EQ.FRST)J=2 IF(I.LT.LCAP2)CA(I,J)=CA(I,J)+CNT IF(I.NE.FRST)THEN J=1 IF(I.LE.FRST+1)J=2 CM(I,J)=CM(I,J)+CNT IF(J.EQ.1.AND.MOD(CAP(I-1),3).GT.0)CMC(I)=CMC(I)+CNT ENDIF ELSE IF(I.LT.LCAP2)THEN J=1 IF(I.LE.FRST+1)J=2 Z(I,J)=Z(I,J)+CNT IF(J.EQ.1.AND.MOD(CAP(I-1),3).GT.0)ZC(I)=ZC(I)+CNT ENDIF ENDIF 130 CONTINUE J=1 IF(LCAP.EQ.FRST)J=2 IF(CAP(LCAP).EQ.1)CN(LCAP,J)=CN(LCAP,J)+CNT C GO TO 40 C *** B-TABLE FORMAT 140 IF(FMT(1:1).NE.'(')THEN DO 180 J=1,NT READ(1,*,ERR=900)(DTABLE(I,J),I=1,NT) 180 IREC=IREC+1 READ(1,*)(NN(I),I=1,NT) READ(1,*)(NS(I),I=1,NT) ELSE DO 181 J=1,NT READ(1,FMT,ERR=900)(DTABLE(I,J),I=1,NT) 181 IREC=IREC+1 READ(1,FMT)(NN(I),I=1,NT) READ(1,FMT)(NS(I),I=1,NT) ENDIF IREC=IREC+2 GO TO 200 C *** TOTALS FORMAT 195 WRITE(7,FMT='(/10X,1Hi,5X,6Hmarked,5X,6Hcaught,5X,8Hreleased, ,5X,11H recaptured,7X,1Hz)') DO 197 I=1,NT IF(FMT(1:1).NE.'(')THEN READ(1,*,ERR=900)NM(I),NN(I),NS(I),R(I),ZZ(I) ELSE READ(1,FMT,ERR=900)NM(I),NN(I),NS(I),R(I),ZZ(I) ENDIF IF(T(I).LE..0D0)T(I)=.1D1 IF(I.GT.1)YR(I)=YR(I-1)+T(I-1) 197 WRITE(7,198)I,NM(I),NN(I),NS(I),R(I),ZZ(I) 198 FORMAT(I11,2F11.0,F13.0,F13.0,F13.0) GO TO 300 C C @@END OF CAPTURE-HISTORY RECORDS PGM LANDS HERE... 200 DO 210 I=1,NT NM(I)=.0D0 DO 205 J=1,I 205 NM(I)=NM(I)+DTABLE(I,J) 210 NU(I)=NN(I)-NM(I) C *** WRITE OUT TABLE WRITE(7,220)IREC 220 FORMAT(/1X,I6,' Input records read'/) ZZ(1)=.0D0 DO 225 I=2,NT IF(IFMT.LE.2)THEN R(I-1)=CA(I-1,1)+CA(I-1,2) ZZ(I)=Z(I,1)+Z(I,2) ELSEIF(IFMT.EQ.3)THEN R(I-1)=.0D0 DO 224 J=I,NT 224 R(I-1)=R(I-1)+DTABLE(J,I-1) ZZ(I)=ZZ(I-1)+R(I-1)-NM(I) ENDIF IF(T(I-1).LE.0.0)T(I-1)=1.0 225 YR(I)=YR(I-1)+T(I-1) IF(YR(NT).NE.NT+YR(1)-1)THEN DO 226 I=1,NT 226 YR(I)=I ENDIF IF(INDEX(MDLS,'AX').LE.0.AND.FLAGAX.EQ.1)WRITE(7,227) 227 FORMAT(/' *** DATA CONTAIN BETWEEN PERIOD RESIGHTING CODES ', / 'WHICH WILL BE IGNORED ***'/ / ' *** TO USE THESE DATA, SELECT MODEL AX ', / ' ***'/) WRITE(7,228)(YR(I),I=1,NT) 228 FORMAT(/' Data summarized in "B-Table" format (See Leslie, ', 1'Chitty and Chitty 1953; Biometrika 40:137-169).'/ 2 /' Time of !'/' last !',10X,'Time of recapture'/ /' capture !',(49I5)) WRITE(7,FMT='(1X,22A5)')(DASH(1),I=1,MIN(NT+2,22)) WRITE(7,250)YR(1),(INT(DTABLE(J,1)),J=1,NT) ZP(NT)=0 DO 230 I=2,NT ZP(NT-I+1)=ZP(NT-I+2)-ZZ(NT-I+2)+ZZ(NT-I+1)+NU(NT-I+2) 230 WRITE(7,250)YR(I),(INT(DTABLE(1,1)),J=1,I-1), , (INT(DTABLE(J,I)),J=I,NT) WRITE(7,*) WRITE(7,FMT='(10H Marked !,49I5)')(INT(NM(J)),J=1,NT) WRITE(7,FMT='(10H Unmarked!,49I5)')(INT(NU(J)),J=1,NT) WRITE(7,FMT='(10H Caught !,49I5)')(INT(NN(J)),J=1,NT) WRITE(7,FMT='(10H Released!,49I5)')(INT(NS(J)),J=1,NT) WRITE(7,FMT='(/ 21H Other summary stats:/)') WRITE(7,FMT='(10H R(i) !,49I5)')(INT(R(J)),J=1,NT) WRITE(7,FMT='(10H z(i) !,49I5)')(INT(ZZ(J)),J=1,NT) WRITE(7,FMT='(10H z''(i) !,49I5)')(INT(ZP(J)),J=1,NT) 250 FORMAT(1X,I7,' !',49I5) 300 READ(TITLE,FMT='(16A8)')TITL RETURN 900 PRINT *,'Error reading record ',IREC,' from input file.' IERR=IERR+1 IF(IERR.LE.10)GO TO 40 STOP END C$CONTROL SEGMENT 'BTBLE' SUBROUTINE GETARY(RCD,X) PARAMETER (IDIM=50) DOUBLE PRECISION X(IDIM) CHARACTER*250 RCD LOGICAL DP I=1 J=1 M=0 X(1)=.0D0 DP=.FALSE. 5 IF(RCD(J:J).EQ.',')THEN IF(M.GT.1)X(I)=X(I)/(10**(M-1)) I=I+1 X(I)=.0D0 DP=.FALSE. M=0 ELSEIF(RCD(J:J).EQ.'.')THEN DP=.TRUE. ELSEIF(RCD(J:J).GE.'0'.AND.RCD(J:J).LE.'9')THEN X(I)=X(I)*10.+ICHAR(RCD(J:J))-ICHAR('0') ELSE GO TO 20 ENDIF IF(DP)M=M+1 J=J+1 IF(J.LE.250)GO TO 5 20 IF(M.GT.1)X(I)=X(I)/(10**(M-1)) RETURN END C$CONTROL SEGMENT 'BTBLE' FUNCTION INUM(STR) CHARACTER*250 STR I=1 INUM=0 5 IF(STR(I:I).LT.'0'.OR.STR(I:I).GT.'9')GO TO 10 INUM=INUM*10+ICHAR(STR(I:I))-ICHAR('0') I=I+1 GO TO 5 10 RETURN END C$CONTROL SEGMENT 'GOF' FUNCTION CHIPRB(IDF,CHSQ) C FUNCTION COMPUTES THE Prob. OF THE Chi-square VALUE 'CHSQ' WITH C 'DF' Degrees of freedom. C C SUBROUTINES NEEDED: C FUNCTION ALGAMA C DOUBLE PRECISION CHIPRB,CHSQ CHIPRB=1. IF(CHSQ .LE. 0.002 .OR. IDF .LT. 1) GO TO 50 A=FLOAT(IDF)/2. X=.5*CHSQ TERM=1. SUM=0. COFN=A IF(X.GE.A.AND.X.GE.13.)GO TO 20 C CONVERGENT SERIES FOR X .LT. A OR .LT. 13. FACT=-A 10 TEMP=SUM SUM=SUM+TERM COFN=COFN+1. TERM=TERM*X/COFN IF(SUM-TEMP)40,40,10 C ASYMPTOTIC SERIES FOR X .GTE. A AND X .GTE. 13. 20 CHIPRB=0. FACT=X 30 TEMP=SUM SUM=SUM+TERM COFN=COFN-1. RATIO=COFN/X TERM=TERM*RATIO IF(SUM.GT.TEMP)GO TO 30 C 40 ARGMNT=ALOG(SUM)-X+A*ALOG(X)-ALGAMA(A) IF(ARGMNT .LT. -100.)GO TO 50 CHIPRB=CHIPRB+EXP(ARGMNT)/FACT 50 RETURN END C$CONTROL SEGMENT 'GOF' FUNCTION ALGAMA(A) C FUNCTION COMPUTES THE LOG OF THE GAMMA FUNCTION GIVEN C THAT THE ARGUMENT IS AN INTEGER OR A MULTIPLE OF 1/2. C ALGAMA=0.0 I=INT(A) B=FLOAT(I) IF(B.NE.A)GO TO 20 IF(A.LE. 2.0)GO TO 40 DO 10 J=2,I-1 10 ALGAMA=ALGAMA+ALOG(FLOAT(J)) GO TO 40 C A IS A MULTIPLE OF 1/2 20 ALGAMA=0.5723649429 IF(A.LE. 0.5)GO TO 40 DO 30 J=1,I*2-1,2 30 ALGAMA=ALGAMA+ALOG(FLOAT(J)/2.) 40 RETURN END C$CONTROL SEGMENT 'GOF' FUNCTION ZPROB(ZVAL) C FUNCTION ZPROB COMPUTES THE PROBABILITY OF A Z VALUE LESS THAN OR C EQUAL TO THE INPUT VALUE FOR A STANDARD NORMAL DISTRIBUTION. C (I.E. P(Z<=Z) ) C THIS PROGRAM WAS CONVERTED FROM THE TEXAS INSTRUMENTS APPLIED C STATISTICS MODULE PROGRAM(#19). C ZPROB=0.0 IF(ABS(ZVAL).GT.5.0)GO TO 10 A=1./(ABS(ZVAL)*.2316419+1.) ZPROB=(A*.31938153-.3565637822*A*A+1.781477937*A*A*A- - 1.821255978*A*A*A*A+1.330274429*A*A*A*A*A)/ / SQRT(EXP(ZVAL*ZVAL)*6.283185307) 10 IF(ZVAL.LT.0.)ZPROB=1.0-ZPROB RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE DEFINE WRITE(7,10)CHAR(12) 10 FORMAT(A1,24X,'P R O G R A M J O L L Y'/ 1 25X,'============= ========= '/ 2 25X,'LAST REVISED : 01/24/91'/) WRITE(7,20) 20 FORMAT(/' Program JOLLY output contains analyses for as many', 1' as 5 different capture-recapture models.'// 2' Model A is the standard Jolly-Seber model for open populations', 3' (Jolly 1965; Biometrika 52: 225-241).'// 460H Model A' is the "Death But No Immigration" model of Darroch, 5' (1959: Biometrika 46: 344-349) and'/ 610X,'Jolly (1965; Biometrika 52: 241-242).'// 7' Model 2 is the capture-resighting model of Brownie and Robson', 8' (1983; Biometrics 39:437-453) permitting an effect'/10X,'of ', 9'initial capture and tagging on first period survival rate.') WRITE(7,30) 30 FORMAT(/' Model B is the Jolly-Seber model with survival', ,' rate assumed constant per unit time'/ /10X,' and time-specific capture probability.'/ //' Model D is the Jolly-Seber model with both survival rate', ,' and capture probability assumed constant per unit time.'//) WRITE(7,40) 40 FORMAT(' ******** Definitions and Notation ********************', 1'***********************************'/) WRITE(7,50) 50 FORMAT(' M(i)',8X,'= Estimated number of marked animals in ', 3 'the population at time i.'// 5 ' p(i)',8X,'= Estimated probability that an animal ', 6'alive at time i is captured in the i-th sample'// 7 ' N(i)',8X,'= Estimated population size at time i'// 8 ' PHI',9X,'= Estimated probability that an animal ', 9'alive at time i survives to time i+1'/) WRITE(7,60) 60 FORMAT(' PHI*',8X,'= Estimated probability that an animal ', 1'caught for the first time in sample i'/ 2 14X,'survives to time i+1 (Model 2 only)'// 3 ' B(i)',8X,'= Estimated number of new animals recruit', 4'ed during the interval i to i+1'/ 5 14X,'and alive at time i+1 (Recruitment includes', 6' birth and immigration)'// 7 ' SE(x)',7X,'= Standard error of parameter x including', 8' non-sampling error terms'/) WRITE(7,70) 70 FORMAT(' SE''(x)',6X,'= Standard error of parameter x ', 1'excluding non-sampling error terms'/ 6 /' COV(X(i,j)) = Covariance between estimates X(i) and', 7' X(j).'/) WRITE(7,90) 90 FORMAT(' r(i)',8X,'= Number of animals caught in sample i', ,', and recaptured later'// /' z(i)',8X,'= Number of animals caught before and after', ,' sample i, but not caught in sample i'// /' z''(i)',7X,'= z(i) + animals caught for the first time ', ,'subsequent to sample i'//) WRITE(7,100) 100 FORMAT(' Notes: Estimates of PHI,N,M,and B under models A and ', ,'A'' are computed using bias corrected formulae'/ /8X,'of Seber (1973:204; The Estimation of Animal ', ,'Abundance and Related Parameters, Griffin, London)'/ /8X,'except when z(i)=0;then N(i)= Number caught in sample i.'/ //8X,'95% Confidence intervals are computed using the ', ,'standard error which includes sampling variation (SE(x)).'/ //8X,'Estimates which cannot be computed due to poor data (m(i),', ,'r(i),z(i),... = zero) are denoted by "Div/0".'/8X,'(Means do', ,' not include these values.)'/) C RETURN END C$CONTROL SEGMENT 'GOF' SUBROUTINE GFIT(Z,R,M,MOD1A,CHISQ,DF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) INTEGER YR,DF,DFLOST DOUBLE PRECISION EX(IDIM,IDIM),M(IDIM),R(IDIM),Z(IDIM), , OBS(IDIM,IDIM),NM,NU,NN,NS EQUIVALENCE (EX,OBS) COMMON /BLK0/NYRS,YR(IDIM),TITLE(16) COMMON /BLK3/BTBLE(IDIM,IDIM),NM(IDIM),NU(IDIM),NN(IDIM),NS(IDIM) C----------------------------------------------------------------------- C PROGRAM TESTS THE DATA FOR Goodness-of-fit TO THE JOLLY MODEL BY C COMPARING THE OBSERVED BTABLE TO THE EXPECTED VALUES. (COLS ARE C POOLED WHEN EXPECTED VALUES ARE SMALL.) C C VARIABLES: M(I)-# MARKED ANIMALS IN THE I-TH SAMPLE. C R(I)-# ANIMALS RELEASED IN TIME I & CAUGHT LATER. C Z(I)-# ANIMALS MARKED before TIME I & CAUGHT LATER. C Z(I)-# Z(I) + SUM(UNMARKED CAUGHT after I). C S(I)-# RELEASED FROM THE I-TH SAMPLE. C C SUBROUTINES NEEDED: C FUNCTION "CHIPRB". C---------------------------------------------------------------------- C INITIALIZE NM1=NYRS-1 CHISQ=.0D0 DFLOST=0 DO 10 I=1,NM1 DO 10 J=I+1,NYRS EX(J,I)=.0D0 10 OBS(I,J)=BTBLE(J,I) WRITE(7,20) 20 FORMAT(//' Observed values (small expected values pooled)'/) C C COMPUTE EXPECTED VALUE MARTIX DO 80 I=1,NM1 IP1=I+1 SUMEXP=.0D0 IF(R(I).LE..0D0)GO TO 35 TMP=R(I)/(Z(I)+R(I)) SUMEXP=M(IP1)*TMP EX(IP1,I)=SUMEXP IP2=I+2 IF(IP2.GT.NYRS)GO TO 80 DO 30 J=IP2,NYRS IF(Z(J-1).LE..0D0)GO TO 35 TMP=TMP*Z(J-1)/(Z(J-1)+R(J-1)) EX(J,I)=TMP*M(J) 30 SUMEXP=SUMEXP+EX(J,I) C ********POOL RECAPTURES & EXPECTED VALUES********* C ** DATA IS POOLED INTO THE 2ND LAST COLUMN ** C ** STARTING AT THE 3RD LAST COL AND PROGRESSING ** C ** TO THE LEFT UNTIL THE EXPECTED VLAUE IN THE ** C ** 2ND LAST COL IS > OR = 2.0. ** C ** REMAINING CELLS WHICH HAVE EXPECTED VALUES ** C ** LESS THAN 2.0 ARE POOLED TO THE ADJACENT ** C ** CELL ON THE LEFT. ** C ************************************************** 35 IF(IP1.GT.NM1)GO TO 80 DO 60 J=IP1,NM1 L=NYRS-J+I IF(EX(NYRS,I).GE..2D1)GO TO 50 EX(NYRS,I)=EX(NYRS,I)+EX(L,I) OBS(I,NYRS)=OBS(I,NYRS)+OBS(I,L) 40 EX(L,I)=.0D0 OBS(I,L)=.0D0 DFLOST=DFLOST+1 GO TO 60 50 IF(EX(L,I).GE..2D1)GO TO 60 LM1=L-1 IF(LM1.LE.I)GO TO 80 EX(LM1,I)=EX(LM1,I)+EX(L,I) OBS(I,LM1)=OBS(I,LM1)+OBS(I,L) GO TO 40 60 CONTINUE C C WRITE MATRICES & COMPUTE Chi-square 80 WRITE(7,90)(OBS(I,J),J=I+1,NYRS) 90 FORMAT(19(1X,F6.1)) WRITE(7,100) 100 FORMAT(//' Expected value matrix (small values pooled)'/) DO 120 I=1,NM1 WRITE(7,90)(EX(J,I),J=I+1,NYRS) DO 110 J=I+1,NYRS IF(EX(J,I).GT..0D0)THEN OBS(I,J)=(OBS(I,J)-EX(J,I))*(OBS(I,J)-EX(J,I))/EX(J,I) CHISQ=CHISQ+OBS(I,J) ELSE OBS(I,J)=.0D0 ENDIF 110 CONTINUE 120 CONTINUE WRITE(7,121) 121 FORMAT(/' Chi-square matrix'/) DO 122 I=1,NM1 122 WRITE(7,90)(OBS(I,J),J=I+1,NYRS) DF=(NYRS-2)*(NYRS-3)/2 - DFLOST IF(MOD1A.EQ.0)GO TO 129 C COMPUTE OBS & EXP OF # UNMARKED IN EACH SAMPLE DO 123 I=2,NYRS OBS(1,I)=NU(I) EX(I,1) =NN(I) DO 123 J=1,I-1 123 EX(I,1)=EX(I,1)*Z(J)/(Z(J)+R(J)) DO 125 I=2,NYRS-1 IF(EX(I,1).LT..2D1)THEN EX(I+1,1)=EX(I+1,1)+EX(I,1) OBS(1,I+1)=OBS(1,I+1)+OBS(1,I) DFLOST=DFLOST+1 ELSE CHISQ=CHISQ+(OBS(1,I)-EX(I,1))**2/EX(I,1) ENDIF 125 CONTINUE WRITE(7,126)(OBS(1,I),I=2,NYRS) 126 FORMAT(/' Unmakred observed:'/(1X,10F12.1)) WRITE(7,127)(EX(I,1),I=2,NYRS) 127 FORMAT(/' Unmarked expected:'/(1X,10F12.1)) DFLOST=(NYRS-1)*(NYRS-2)/2 - DFLOST C PRINT TOTAL Chi-square, DF, AND Prob. 129 WRITE(7,130)CHISQ,DF 130 FORMAT(//40X,40('-')/ 1 40X,'! Chi-square value ',F10.4,10X,'!'/ 2 40X,'! Degrees of freedom ',I8,10X,'!') IF(DF.GT. 0)THEN PROB=CHIPRB(DF,CHISQ) WRITE(7,140)PROB 140 FORMAT( 40X,'! Prob. of a larger chi-square ',F7.4,' !'/ 1 40X,40('-')/) ELSE WRITE(7,160) 160 FORMAT(40X,'! INSUFFICIENT DATA FOR GOF TEST',7X,'!'/40X,40('-')/) DF=0 ENDIF RETURN END C$CONTROL SEGMENT 'GOF' SUBROUTINE GOF2(M1A,ZP,U,TSUM2,I2,IDFA,CHIA) C C THIS SUBROUTINE PERFORMS A Goodness-of-fit TEST USING C CONTINGENCY TABLES. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) DIMENSION ZP(IDIM),U(IDIM),EX1(2,4),EX2(2,4),OBS1(2,4),OBS2(2,4) INTEGER YR COMMON /BLK0/LAST,YR(IDIM),TITLE(16) COMMON /BLK2/CA(IDIM,2),CM(IDIM,2),CN(IDIM,2),Z(IDIM,2),CMC(IDIM), , ZC(IDIM) TSUM1=0. TSUM2=0. A1M=1-M1A I1=0 I2=0 IDFA=0 WRITE(7,10) 10 FORMAT(/41X,'Contingency table Goodness-of-fit test'/ / 41X,'--------------------------------------'/ / /20X,'Component 1',55X,'Component 2'/) IF(CM(2,2).LE..0D0)GO TO 120 C C FOR EACH SAMPLE PERIOD FROM 2,3,...LAST-1 C COMPUTE 2 EXPECTED VALUE MATRICES AND C Chi-square VALUES. ALSO, PRINT OUT CONTINGENCY C TABLES C DO 90 I=2,LAST-1 OBS2(1,1)=CA(I,1) OBS2(1,2)=CA(I,2) OBS2(2,1)=CN(I,1) OBS2(2,2)=CN(I,2) OBS2(1,3)=.0D0 OBS2(2,3)=.0D0 OBS2(1,4)=.0D0 OBS2(2,4)=.0D0 OBS1(1,1)=CM(I,1)-CMC(I) OBS1(1,2)=CMC(I) OBS1(1,3)=CM(I,2) OBS1(1,4)=M1A*U(I) OBS1(2,1)=Z(I,1)-ZC(I) OBS1(2,2)=ZC(I) OBS1(2,3)=Z(I,2) OBS1(2,4)=M1A*(ZP(I)-Z(I,1)-Z(I,2)) CALL EXVAL(OBS2,EX2) CALL EXVAL(OBS1,EX1) C IF(M1A.LE.0)WRITE(7,20)YR(I) 20 FORMAT(/88X,'1st cap before i-1,!'/' I=',I4,25X,2(' 1st cap'), ,40X,'not cap cap ! 1st cap'/ /33X,'before i in i',40X,2(' in i-1'),' ! in i-1') IF(M1A.NE.0)WRITE(7,30)YR(I) 30 FORMAT(/88X,'1st cap before i-1,!'/ /' i=',I4,27X,'1st cap 1st cap', ,40X,'not cap cap in ! 1st cap 1st cap'/ /33X,'before i in i',42X,'in i-1 in i-1 ! in i-1 after i-1') IF(M1A.LE.0)WRITE(7,40) IF(M1A.GT.0)WRITE(7,41) 40 FORMAT(32X,'+',2('--------+'),38X,'+',3('--------+')) 41 FORMAT(32X,'+',2('--------+'),38X,'+',4('--------+')) WRITE(7,50)(OBS2(1,J),J=1,2),(OBS1(1,J),J=1,3+M1A) 50 FORMAT(' Cap in i, released & recap !',2(F8.2,'!'), , 15X,' Cap in i !',4(F8.2,'!')) WRITE(7,60)(EX2(1,J),J=1,2),(EX1(1,J),J=1,3+M1A) 60 FORMAT(' Expected value',17X,'!',2(F8.2,'!'),15X, , ' Expected value',8X,'!',4(F8.2,'!')) IF(M1A.LE.0)WRITE(7,40) IF(M1A.GT.0)WRITE(7,41) WRITE(7,70)(OBS2(2,J),J=1,2),(OBS1(2,J),J=1,3+M1A) 70 FORMAT(' Cap in i, released & not recap !',2(F8.2,'!'), , 15X,' Cap after, not in i !',4(F8.2,'!')) WRITE(7,60)(EX2(2,J),J=1,2),(EX1(2,J),J=1,3+M1A) C IDF1=0 SUM1=.0D0 IDF2=0 SUM2=.0D0 IF(EX2(1,1).GE..2D1.AND.EX2(1,2).GE..2D1.AND. . EX2(2,1).GE..2D1.AND.EX2(2,2).GE..2D1)THEN IDF2=IDF2+1 DO 71 II=1,2 DO 71 J=1,2 71 SUM2=SUM2+(OBS2(II,J)-EX2(II,J))**2/EX2(II,J) ENDIF C POOL TABLE ( BY COLS) UNTIL ALL EXP. VALS >= 2.0 CALL POOL(EX1,OBS1) DO 73 J=1,4 IF(EX1(1,J).LT..2D1.OR.EX1(2,J).LT..2D1)GO TO 73 IDF1=IDF1+1 DO 72 II=1,2 72 SUM1=SUM1+(OBS1(II,J)-EX1(II,J))**2/EX1(II,J) 73 CONTINUE IDF1=MAX(IDF1-1,0) IF(IDF1.LT.1)SUM1=0. C C OUTPUT TOTAL Chi-square OF EACH TABLE & Prob. C P1=CHIPRB(IDF1,SUM1) P2=CHIPRB(IDF2,SUM2) IF(M1A.LE.0)WRITE(7,40) IF(M1A.GT.0)WRITE(7,41) WRITE(7,80)SUM2,IDF2,P2,SUM1,IDF1,P1 80 FORMAT(2(' Chi-square = ',F8.4,' with ',I2, , ' Degrees of freedom (Prob.=',F6.4,') ')) TSUM1=TSUM1+SUM1 TSUM2=TSUM2+SUM2 I1=I1+IDF1 I2=I2+IDF2 90 CONTINUE WRITE(7,100) 100 FORMAT(/1X,132('-')) P1=CHIPRB(I1,TSUM1) P2=CHIPRB(I2,TSUM2) WRITE(7,80)TSUM2,I2,P2,TSUM1,I1,P1 CHIA=TSUM1+TSUM2 IDFA=I1+I2 PROB=CHIPRB(IDFA,CHIA) WRITE(7,110)CHIA,IDFA,PROB 110 FORMAT(25X,'Overall chi-square = ',F8.4,' with ',I2, ,' degrees of freedom (Prob. =',F6.4,')'/1X,132('-')) GO TO 140 120 WRITE(7,130) 130 FORMAT(/' **** INSUFFICIENT DATA FOR THIS TEST ****'/) 140 RETURN END C$CONTROL SEGMENT 'GOF' SUBROUTINE EXVAL(OBS,E) C C THIS SUBROUTINE COMPUTES THE Chi-square VALUE FOR C A 2 X 3 CONTINGENCY TABLE. C CONTINGENCY TABLE IS DESCRIBED IN SUBROUTINE "GOF2" C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION OBS(2,4),OBS1(2,4),E(2,4),ROW(2),COL(4) ROW(1)=OBS(1,1)+OBS(1,2)+OBS(1,3)+OBS(1,4) ROW(2)=OBS(2,1)+OBS(2,2)+OBS(2,3)+OBS(2,4) TOTL=ROW(1)+ROW(2) IF(TOTL.LE..0D0)GO TO 70 DO 10 J=1,4 OBS1(1,J)=OBS(1,J) OBS1(2,J)=OBS(2,J) COL(J)=OBS(1,J)+OBS(2,J) E(1,J)=ROW(1)*COL(J)/TOTL 10 E(2,J)=ROW(2)*COL(J)/TOTL 70 RETURN END C$CONTROL SEGMENT 'GOF' SUBROUTINE POOL(E,OBS) C POOLS 2X3 OR 2X4 CONTINGENCY TABLE C 9 DIFFERENT WAYS OF POOLING ARE POSSIBLE FROM A 2X4 TABLE. EACH C OF THESE WAYS IS TRIED TO SEE IF THERE IS A WAY TO ONLY REDUCE C THE NUMBER OF COLUMNS BY 1. C DOUBLE PRECISION E(2,4),OBS(2,4),EE(2,5),EMIN INTEGER II(9),JJ(9),KK(9) DATA II/1,1,1,2,2,3,1,1,2/,JJ/2,3,4,3,4,4,2,2,3/,KK/6*5,3,4,4/ DATA EE/10*.0D0/ C FIRST, SEE IF ANY POOLING IS NECESSARY IF(EMIN(E).GT..2D1)RETURN C NOW, FOR EACH WAY OF POOLING, POOL & CHECK IF ALL EXPECTED C VALUES ARE > 2. DO 10 IWAY=1,9 C MAKE A WORKING COPY OF E. DO 5 I=1,4 EE(1,I)=E(1,I) 5 EE(2,I)=E(2,I) I=II(IWAY) J=JJ(IWAY) K=KK(IWAY) IF(J.EQ.4.OR.IWAY.GT.4)GO TO 10 C POOL COLS I & J (& K IF IWAY>6) EE(1,I)=EE(1,I)+EE(1,J)+EE(1,K) EE(2,I)=EE(2,I)+EE(2,J)+EE(2,K) EE(1,J)=.0D0 EE(1,K)=.0D0 EE(2,J)=.0D0 EE(2,K)=.0D0 IF(EMIN(EE).GE..2D1)GO TO 20 10 CONTINUE C AT THIS POINT, POOLING DOESN'T HELP...SO RETURN RETURN C POOL EXP & OBS TABLES 20 E(1,I)=E(1,I)+E(1,J) E(2,I)=E(2,I)+E(2,J) OBS(1,I)=OBS(1,I)+OBS(1,J) OBS(2,I)=OBS(2,I)+OBS(2,J) E(1,J)=.0D0 E(2,J)=.0D0 OBS(1,J)=.0D0 OBS(2,J)=.0D0 IF(IWAY.GT.6)THEN E(1,I)=E(1,I)+E(1,K) E(2,I)=E(2,I)+E(2,K) OBS(1,I)=OBS(1,I)+OBS(1,K) OBS(2,I)=OBS(2,I)+OBS(2,K) E(1,K)=.0D0 E(2,K)=.0D0 OBS(1,K)=.0D0 OBS(2,K)=.0D0 ENDIF RETURN END C$CONTROL SEGMENT 'GOF' FUNCTION EMIN(EE) C FINDS THE MINIMUN NON-ZERO VALUE IN A 2X4 CONTINGENCY TABLE DOUBLE PRECISION EE(2,5),EMIN EMIN=.1D55 DO 5 I=1,4 IF(EE(1,I).LT.EMIN.AND.EE(1,I).GT..0D0)EMIN=EE(1,I) IF(EE(2,I).LT.EMIN.AND.EE(2,I).GT..0D0)EMIN=EE(2,I) 5 CONTINUE RETURN END C$CONTROL SEGMENT 'GOF' SUBROUTINE GOFM2(Z,M,TOTCHI,ITOTDF) C C SUBROUTINE TO COMPUTE A Chi-square VALUE WHICH TESTS C THE HYPOTHESIS THAT THE DATA FIT MODEL 2. C NOTATION: LOWER CASE V(I,J) = RE-RECAPTURE DATA C = BTABLE(J,I) C LOWER CASE U(I,J) = N(I,J)-V(I,J) C = BTABLE(I,J)-BTABLE(J,I) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) INTEGER YR DIMENSION MSTAR(IDIM,IDIM),EX(IDIM,3),OBS(IDIM,3),Z(IDIM) DOUBLE PRECISION M(IDIM),V(IDIM),U(IDIM) COMMON /BLK0/K,YR(IDIM),TITLE(16) COMMON /BLK2/CA(IDIM,2),CM(IDIM,2),CN(IDIM,2),ZZ(IDIM,2), , CMC(IDIM),ZC(IDIM) COMMON /BLK3/BTABLE(IDIM,IDIM) EQUIVALENCE (V(1),CA(1,1)),(U(1),CA(1,2)) TOTCHI=0. ITOTDF=0 WRITE(7,10)TITLE 10 FORMAT(/1X,16A8//45X,'Goodness-of-fit test for model 2'/ /45X,32('-')/' Contingency tables of pooled observed & expected', ,' values as specified in Brownie & Robson paper.') C COMPUTE MSTAR(I,J) NEEDED FOR MODEL 2 GOF TEST MSTAR(2,1)=BTABLE(2,1) DO 20 I=3,K MSTAR(I,1)=BTABLE(I,1) DO 20 J=2,I-1 20 MSTAR(I,J)=MSTAR(I,J-1)+BTABLE(I,J) C C COMPUTE CONTINGENCY TABLES & Chi-square VALUES C DO 110 I=2,K-2 IF(Z(I).LT..2D1.OR.V(I).LT..2D1.OR.U(I).LT..2D1)GO TO 110 C COMPUTE TABLE OF OBSERVED & EXPECTED VALUES DO 30 J=1,K-I IPJ=I+J OBS(J,1)=MSTAR(IPJ,I-1) OBS(J,2)=BTABLE(I,IPJ) OBS(J,3)=BTABLE(IPJ,I)-BTABLE(I,IPJ) TMP=(MSTAR(IPJ,I))/(Z(I+1)+M(I+1)) EX(J,1)=Z(I)*TMP EX(J,2)=V(I)*TMP 30 EX(J,3)=U(I)*TMP C C POOL ROWS IF EXPECTED VALUE < 2.0 C POOL TO THE ROW WITH THE SMALLEST EXPECTED C VALUE FOR EACH COLUMN. C DO 40 J=1,K-I-1 IF(FINDM(K-I,EX,-1,NPTR).GE..2D1)GO TO 50 TMP=FINDM(K-I,EX,NPTR,N2PTR) IF(N2PTR*NPTR.GT.0)CALL POOLR(EX,OBS,NPTR,N2PTR) 40 CONTINUE C C COMPUTE Chi-square FOR CONTINGENCY TABLE 50 CHI=0. IDF=-2 WRITE(7,60)YR(I) 60 FORMAT(/' i=',I4,1X,'j',4X,'m*(i+j,i-1)',7X,'v(i,j)',6X, , 'u(i,j)',20X,'Expected values') DO 80 J=1,K-I IF(EX(J,1).EQ..0D0)GO TO 80 DO 70 L=1,3 70 CHI=CHI+(EX(J,L)-OBS(J,L))*(EX(J,L)-OBS(J,L))/EX(J,L) IDF=IDF+2 WRITE(7,90)YR(I+J),(OBS(J,L),L=1,3),(EX(J,L),L=1,3) 80 CONTINUE 90 FORMAT(7X,I4,2(3X,F10.4,5X,F10.4,2X,F10.4,4X)) C IF(IDF.LE.0)GO TO 110 PROB=CHIPRB(IDF,CHI) WRITE(7,100)CHI,IDF,PROB 100 FORMAT(/' Chi-square = ',F10.4,' WITH ',I4,' degrees of', , ' freedom : Prob = ',F10.4) TOTCHI=TOTCHI+CHI ITOTDF=ITOTDF+IDF 110 CONTINUE C C COMPUTE & PRINT PROB OF TOTAL Chi-square C PROB=CHIPRB(ITOTDF,TOTCHI) WRITE(7,120)TOTCHI,ITOTDF,PROB 120 FORMAT(/20X,'Total chi-square = ',F12.4/ / 20X,'Degrees of freedom = ',I7/ / 20X,'Probability = ',F12.4) C RETURN END C$CONTROL SEGMENT 'GOF' FUNCTION FINDM(NPER,ARY,ISELF,IROW) C C SUBPROGRAM TO FIND THE MINIMUM EXPECTED VALUE FOR C A SPECIFIC COLUMN IN THE CONTINGCY TABLE (OTHER THAN C IT SELF). IROW POINTS TO THE ROW IN THE TABLE WHERE C THE MINIMUM VALUE IS. C PARAMETER (IDIM=50) DOUBLE PRECISION FINDM,X,ARY(IDIM,3) FINDM=.9D66 IROW=0 DO 10 J=1,3 DO 10 I=1,NPER X=ARY(I,J) IF(I.EQ.ISELF.OR.X.LE..0D0.OR.X.GE.FINDM)GO TO 10 FINDM=X IROW=I 10 CONTINUE RETURN END C$CONTROL SEGMENT 'GOF' SUBROUTINE POOLR(EX,OBS,ROW1,ROW2) C C SUBPROGRAM POOLS EXPECTED & OBSERVED VALUES FROM ROW1 C OF CONTINGENCY TABLE INTO ROW2. C PARAMETER (IDIM=50) INTEGER ROW1,ROW2 DOUBLE PRECISION EX(IDIM,3),OBS(IDIM,3) C DO 10 I=1,3 EX(ROW2,I)=EX(ROW2,I)+EX(ROW1,I) OBS(ROW2,I)=OBS(ROW2,I)+OBS(ROW1,I) EX(ROW1,I)=.0D0 10 OBS(ROW1,I)=.0D0 RETURN END C$CONTROL SEGMENT 'M1A' SUBROUTINE M1A(T,NU,NN,NS,R,ZPRIME,MDL,CHIA,IDFA,PRT) C C SUBROUTINE TO COMPUTE ESTIMATES FROM JOLLY MODEL WITH C DEATH AND NO IMMIGRATION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) DOUBLE PRECISION T(IDIM),NU(IDIM),NN(IDIM),NS(IDIM),R(IDIM), , ZPRIME(IDIM) INTEGER YR,PRT LOGICAL AV1 CHARACTER*6 UNDEF DIMENSION YN(IDIM),YNSD(IDIM),YNSD2(IDIM) COMMON /BLK0/NTIME,YR(IDIM),TITLE(16) DATA YNAVG,YNSDAV,YNSD2A,YPAVG,YFSD2A,YPHIAV,YFISDA/7*.0D0/ DATA AVGT,AVST,YNSD,YNSD2/2*.0D0,IDIM*.1D30,IDIM*.1D30/ DATA UNDEF/' Div/0'/ C NM1=NTIME-1 NM2=NTIME-2 WRITE(7,10)CHAR(12),TITLE 10 FORMAT(A1,16A8) WRITE(7,20) 20 FORMAT(/' -- Model A'' -- Estimates from capture-recapture', 1 ' data with death but no immigration.'/) WRITE(7,30) 30 FORMAT(/8X,'==== Survival rate estimates between sampling ', ,'periods ==== ', ,'Interval Survival rate estimates per unit time (PHI**(1/T(i))'/ /' Period PHI SE(PHI) ',8HSE'(PHI),' 95% Conf. ', ,'Interval COV(PHI(i,i-1)) Length PHI SE(PHI)', ,' 95% Conf. interval COV(PHI(i,i-1))'/1X,132('-')) YN(1)=(NS(1)+.1D1)/(R(1)+.1D1)*ZPRIME(1)+NN(1) YNAVG=YN(1) NPHI=0 NSE=0 AV1=.TRUE. DO 80 I=1,NM2 YN(I+1)=(NS(I+1)+.1D1)/(R(I+1)+.1D1)*ZPRIME(I+1)+NN(I+1) YNAVG=YNAVG+YN(I+1) IF(R(I)*NS(I).GT..0D0)THEN YPHI=YN(I+1)/(NS(I)/R(I)*ZPRIME(I)+NS(I)) NPHI=NPHI+1 YPHIAV=YPHIAV+YPHI PHIT=YPHI**(.1D1/T(I)) IF(YPHI.LE..0D0)PHIT=.0D0 IF(R(I+1)*NS(I+1).GT..0D0)THEN C COMPUTE VARIANCE OF PHI C=.1D1/R(I)-.1D1/NS(I) G=(1.0/R(I+1))-(1.0/NS(I+1)) YE=YN(I+1)-NN(I+1) YVFI2=YE*(YE+NS(I+1))/YN(I+1)/YN(I+1)*G+ 1 (YN(I)-NN(I))/(YN(I)-NN(I)+NS(I))*C YVARFI=(YVFI2+(.1D1-YPHI)/YN(I+1))*YPHI*YPHI NSE=NSE+1 YFISD=DSQRT(YVARFI) YFISDA=YFISDA+YVARFI YVFI2=YVFI2*YPHI*YPHI YFISD2=DSQRT(YVFI2) YFSD2A=YFSD2A+YVFI2 CI1=YPHI-.196D1*YFISD CI2=YPHI+.196D1*YFISD TT=(.1D1-T(I))/T(I) SEPHIT=YFISD*(YPHI**TT/T(I)) AVGT=AVGT+PHIT AVST=AVST+SEPHIT*SEPHIT IF(I.EQ.1)THEN WRITE(7,41)YR(I),YPHI,YFISD,YFISD2,CI1,CI2,T(I),PHIT, , SEPHIT,PHIT-.196D1*SEPHIT,PHIT+.196D1*SEPHIT ELSE YPHIM1=YN(I)/(YN(I-1)-NN(I-1)+NS(I-1)) YFICOV=-YPHIM1*YPHI*(YN(I)-NN(I))/YN(I)*C COVPHT=YFICOV*YPHIM1**((.1D1-T(I-1))/T(I-1))* * YPHI**TT /T(I-1)/T(I) WRITE(7,40)YR(I),YPHI,YFISD,YFISD2,CI1,CI2,YFICOV,T(I), , PHIT,SEPHIT,PHIT-.196D1*SEPHIT,PHIT+.196D1*SEPHIT,COVPHT AVST=AVST+COVPHT*.2D1 YFISDA=YFISDA+YFICOV*.2D1 YFSD2A=YFSD2A+YFICOV*.2D1 ENDIF ELSE IF(I.EQ.1)WRITE(7,42)YR(I),YPHI,(UNDEF,J=1,4),' ',T(I), , PHIT,(UNDEF,J=1,3) IF(I.NE.1)WRITE(7,42)YR(I),YPHI,(UNDEF,J=1,5),T(I), , PHIT,(UNDEF,J=1,4) AV1=.FALSE. ENDIF ELSE IF(I.EQ.1)WRITE(7,43)YR(I),(UNDEF,J=1,5),' ',T(I), , (UNDEF,J=1,4) IF(I.NE.1)WRITE(7,43)YR(I),(UNDEF,J=1,6),T(I), , (UNDEF,J=1,5) AV1=.FALSE. ENDIF IF(MOD(YR(I),5).EQ.0)WRITE(7,40) 80 CONTINUE 40 FORMAT(1X,I4,F7.4,1X,F6.4,F8.4,F11.4,' -',F7.4,F17.10, , 6X,F7.4,2F7.4,F9.4,' -',F7.4,F19.10) 41 FORMAT(1X,I4,F7.4,1X,F6.4,F8.4,F11.4,' -',F7.4,17X, , 6X,F7.4,2F7.4,F9.4,' -',F7.4,F19.10) 42 FORMAT(1X,I4,F7.4,1X,A6,A8,A11,' -',A7,A17, , 6X,F7.4,2(1X,A6),3X,A6,' - ',A6,7X,A12) 43 FORMAT(1X,I4,A7,1X,A6,A8,A11,' -',A7,A17, , 6X,F7.4,2(1X,A6),3X,A6,' - ',A6,7X,A12) IF(NPHI.GT.0)THEN YPHIAV=YPHIAV/(NPHI) AVGT=AVGT/(NPHI) IF(AV1)THEN YFSD2A=DSQRT(YFSD2A)/(NSE) YFISDA=DSQRT(YFISDA)/(NSE) CI1=YPHIAV-.196D1*YFISDA CI2=YPHIAV+.196D1*YFISDA AVST=DSQRT(AVST)/(NSE) WRITE(7,100)YPHIAV,YFISDA,YFSD2A,CI1,CI2, , AVGT,AVST,AVGT-.196D1*AVST,AVGT+.196D1*AVST ELSE WRITE(7,101)YPHIAV,(UNDEF,J=1,4),AVGT,(UNDEF,J=1,3) ENDIF ELSE WRITE(7,102)(UNDEF,J=1,9) ENDIF 100 FORMAT(1X,132('-')/' MEAN ',2(F6.4,1X),1X,F6.4,5X,F6.4,' - ', ,F6.4,30X,2(1X,F6.4),3X,F6.4,' - ',F6.4/1X,132('-')) 101 FORMAT(1X,132('-')/' MEAN ',F6.4,A7,A8,A11,' - ', ,A6,30X,F7.4,A7,A9,' - ',A6/1X,132('-')) 102 FORMAT(1X,132('-')/' MEAN ',A6,A7,A8,A11,' - ', ,A6,30X,A7,A7,A9,' - ',A6/1X,132('-')) WRITE(7,880) 880 FORMAT(//' Period p SE(p)',5X, ,'95% Conf. interval'/1X,52('-')) PMEAN=.0D0 PMVAR=.0D0 NP=0 DO 881 I=2,NM1 IF(R(I)*NN(I).GT..0D0)THEN NP=NP+1 PI=NN(I)/(NN(I)+NS(I)*ZPRIME(I)/R(I)) PMEAN=PMEAN+PI SEP=DSQRT(PI*PI*(.1D1-PI)*(.1D1-PI)* * (.1D1/R(I)-.1D1/NS(I)+.1D1/NN(I)+.1D1/ZPRIME(I))) PMVAR=PMVAR+SEP*SEP CI1=.196D1*SEP CI2=PI+CI1 CI1=PI-CI1 WRITE(7,882)YR(I),PI,SEP,CI1,CI2 882 FORMAT(1X,I4,F12.4,F12.4,F12.4,' - ',F7.4) ELSE WRITE(7,885)YR(I),(UNDEF,J=1,4) 885 FORMAT(1X,I4,A12,A12,A12,' - ',A7) ENDIF 881 IF(MOD(YR(I),5).EQ.0)WRITE(7,882) PMEAN=PMEAN/(NP) PMVAR=DSQRT(PMVAR)/(NP) CI1=.196D1*PMVAR CI2=PMEAN+CI1 CI1=PMEAN-CI1 WRITE(7,884)PMEAN,PMVAR,CI1,CI2 884 FORMAT(1X,52('-')/' MEAN',F12.4,F12.4,F12.4,' - ',F7.4/ /1X,52('-')) WRITE(7,110) 110 FORMAT(//' Period',4X,4H N ,7X,5HSE(N),8X,6HSE'(N), ,10X,'95% confidence interval'/1X,77('-')) 120 FORMAT(1X,I4,F9.2,1X,F12.2,2(1X,F13.2),' - ',F14.2,5X,F14.2) 121 FORMAT(1X,I4,F9.2,A13,2A14,' - ',A14,5X,A14) AV1=.TRUE. DO 130 I=1,NM1 IF(R(I)*NS(I).GT..0D0)THEN DO 60 J=1,I-1 60 YNSDAV=YNSDAV+.2D1*(YN(I)-YN(I)*YN(J)/YN(1)) YNVAR=(YN(I)-NN(I))*(YN(I)-NN(I)+NS(I))* * (.1D1/R(I)-.1D1/NS(I)) YNSD2(I)=DSQRT(YNVAR) IF(YN(I).LT.YN(1))YNVAR=YNVAR+YN(I)-YN(I)*YN(I)/YN(1) YNSD(I)=DSQRT(YNVAR) YNSD2A=YNSD2A+YNSD2(I)*YNSD2(I) YNSDAV=YNSDAV+YNSD(I)*YNSD(I) CI1=YN(I)-.196D1*YNSD(I) CI2=YN(I)+.196D1*YNSD(I) WRITE(7,120)YR(I),YN(I),YNSD(I),YNSD2(I),CI1,CI2 ELSE WRITE(7,121)YR(I),YN(I),(UNDEF,J=1,4) AV1=.FALSE. ENDIF 130 IF(MOD(YR(I),5).EQ.0)WRITE(7,120) YNAVG=YNAVG/(NM1) IF(AV1)THEN YNSD2A=DSQRT(YNSD2A)/(NM1) YNSDAV=DSQRT(YNSDAV)/(NM1) CI1=YNAVG-.196D1*YNSDAV CI2=YNAVG+.196D1*YNSDAV WRITE(7,140)YNAVG,YNSDAV,YNSD2A,CI1,CI2 ELSE WRITE(7,141)YNAVG,(UNDEF,J=1,4) ENDIF 140 FORMAT(1X,73('-')/' MEAN ',F8.2,1X,F12.2,2(1X,F13.2),' - ', ,F14.2/1X,73('-')) 141 FORMAT(1X,73('-')/' MEAN ',F8.2,A13,2A14,' - ',A14/1X,73('-')) IF(PRT.EQ.1)THEN WRITE(7,10)CHAR(10),TITLE WRITE(7,20) WRITE(7,150) 150 FORMAT(/' Covariance(N(i),N(j)) (i=2,s-1 ; j=1,i-1) : '/) DO 160 I=2,NTIME-1 160 WRITE(7,170)YR(I),((YN(I)-YN(I)*YN(J)/YN(1)),J=1,I-1) ENDIF 170 FORMAT(/' ',I4,10F12.2/(5X,10F12.2)) CALL GOF2(1,ZPRIME,NU,TSTAT,ITMP,IDFA,CHIA) IF(IDFA.GT.0)GO TO 200 WRITE(7,190)CHAR(12) 190 FORMAT(A1,' B-TABLE goodness-of-fit test for model A',1H'/) CALL GFIT(ZPRIME,R,NN,MDL,CHIA,IDFA) 200 RETURN END C$CONTROL SEGMENT 'MOD2' SUBROUTINE MODEL2(T,NM,NN,NS,Z,CHI2,IDF2,CC) C C SUBPROGRAM TO COMPUTE ESTIMATES OF CAPTURE Prob. (P), C SURVIVAL RATES OF UNMARKED (S*), AND SURVIVAL RATES C OF MARKED (S) ACCORDING TO BROWNIE & ROBSONS C GENERALIZATION OF THE JOLLY-SEBER MODEL. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) INTEGER YR,CC LOGICAL AV1 CHARACTER*6 UNDEF DOUBLE PRECISION NM(IDIM),NN(IDIM),NS(IDIM),Z(IDIM),S(IDIM), , P(IDIM),T(IDIM),SEP(IDIM),SES(IDIM),SCOV(IDIM) COMMON /BLK0/NTIME,YR(IDIM),TITLE(16) COMMON /BLK2/CA(IDIM,2),CM(IDIM,2),CN(IDIM,2),ZZ(IDIM,2), , CMC(IDIM),ZC(IDIM) DATA AVGS,AVGP,AVGSP,AVGSES,AVGSEP,AVSESP,AVGSN,AVGPN,AVGSPN, , AVSESN,AVSSPN/11*.0D0/,AV1/.TRUE./ DATA AVGT,AVST,S,SES,P,SEP,SCOV/2*.0D0,IDIM*.1D30,IDIM*.1D30, , IDIM*.1D30,IDIM*.1D30,IDIM*.1D30/,UNDEF/' Div/0'/ C WRITE(7,10)CHAR(12),TITLE 10 FORMAT(A1,16A8) WRITE(7,20) 20 FORMAT(/' -- Model 2 -- Brownie & Robsons generalization of', ,' the Jolly-Seber model (model A) with'/ /15X,'survival rates of newly captured animals different from ', ,'survival rates'/15X,'of previously captured animals.') DO 30 I=1,NTIME-2 IF(NN(I).NE.NS(I))THEN WRITE(7,35) RETURN ENDIF 30 CONTINUE 35 FORMAT(/' *** MODEL 2 NOT APPROPRIATE...***'/ / ' *** # CAPTURED > # RELEASED ***'/) CC=1 WRITE(7,40) 40 FORMAT(/' Some additional definitions and notation for', 1' model 2:'/ 2/5X,'v(i,j) = number of animals caught in period j and', 3 'next recaptured in period i.'/) WRITE(7,50) 50 FORMAT(5X,'u(i,j) = Number of animals marked in peroiod j ', 1'and first caught in period i.'/ 2/' m*(i,j) = (u(i,1)+u(i,2)+...+u(i,j))+(v(i,2)+v(i,3)+', 3'...+v(i,j))'/) WRITE(7,60) 60 FORMAT(5X,'u(i) = Number of animals marked in period i and', 1' ecaptured laterR'/ 2/5X,'v(i) = Number of animals caught in period i and ', 3'recaptured later.'/ 4/5X,'m''(i) = Number of marked animals caught in period i'/ 5/5X,'n''(i) = Number of animals captured, marked and ', 6'released in period i.'/ 7/5X,'zz(i) = Number of animals marked before period i that', 8' are not caught in i, but are caught after i.'/) WRITE(7,70) 70 FORMAT(/6X,'Period',8X,'u(i)',10X,'v(i)',10X,'m',1H', 1'(i)',10X,'n',1H','(i)',10X,'zz(i)'/) DO 80 I=1,NTIME WRITE(7,90)YR(I),CA(I,2),CA(I,1),NM(I),NN(I)-NM(I),Z(I) 80 IF(MOD(YR(I),5).EQ.0)WRITE(7,90) 90 FORMAT(8X,I4,F12.0,F14.0,3F15.0) WRITE(7,100) 100 FORMAT(/8X,'== Surv rate ests between sampling ', ,'periods ==',17X, ,'Interval Surv rate ests per unit time (PHI**(1/T(i))'/ /' Period PHI* SE(PHI*) 95% Conf. ', ,'interval',20X,'Length PHI* SE(PHI*)',9X, ,'95% Conf. interval'/1X,132('-')) C C COMPUTE ESTIMATES USING FORMULAE IN BROWNIE C AND ROBSONS PAPER. C S(1)=.0D0 DO 150 I=1,NTIME-2 TN=NN(I)-NM(I) TM=NM(I) TZ=Z(I) TV=CA(I,1) TU=CA(I,2) TM1=NM(I+1) TZ1=Z(I+1) TV1=CA(I+1,1) IF(TV1*TN*(TM1+TZ1).GT..0D0)THEN SP=TM1/TV1*(TZ1+TV1)/(TM1+TZ1)*TU/TN AVGSP=AVGSP+SP AVGSPN=AVGSPN+.1D1 PHIT=.0D0 IF(SP.GT..0D0)PHIT=SP**(.1D1/T(I)) AVGT=AVGT+PHIT IF(TM1*TU.GT..0D0)THEN TEMP=.1D1/TV1-.1D1/(TM1+TZ1)+(TV1-TM1-TZ1)/(TM1*(TV1+TZ1)) IF(DABS(TEMP).LT..1D-8)TEMP=.0D0 SESP=DSQRT(SP*SP*(.1D1/TU-.1D1/TN+TEMP)) AVSESP=AVSESP+SESP*SESP AVSSPN=AVSSPN+.1D1 TT=(.1D1-T(I))/T(I) SEPHIT=.0D0 IF(SP.GT..0D0)SEPHIT=SESP*(SP**TT/T(I)) AVST=AVST+SEPHIT*SEPHIT CI1=SP-.196D1*SESP CI2=SP+.196D1*SESP C WRITE(7,140)YR(I),SP,SESP,CI1,CI2,T(I),PHIT,SEPHIT, , PHIT-.196D1*SEPHIT,PHIT+.196D1*SEPHIT ELSE WRITE(7,141)YR(I),SP,(UNDEF,J=1,3),T(I),PHIT,(UNDEF,J=1,3) AV1=.FALSE. ENDIF ELSE WRITE(7,142)YR(I),(UNDEF,J=1,4),T(I),(UNDEF,J=1,4) AV1=.FALSE. ENDIF IF(MOD(YR(I),5).EQ.0)WRITE(7,140) IF(TV1*TM*(TM1+TZ1).GT..0D0)THEN S(I)=TM1/TV1*(TZ1+TV1)/(TM1+TZ1)*TV/TM AVGSN=AVGSN+.1D1 AVGS=AVGS+S(I) IF(TV.GT..0D0)THEN VARS=S(I)*S(I)*(.1D1/TV-.1D1/TM+TEMP) SES(I)=DSQRT(VARS) AVGSES=AVGSES+VARS IF(I.GT.2)THEN SCOV(I)=-S(I-1)*S(I)*(.1D1/TV-(TM+TZ)/(TM*(TV+TZ))) AVGSES=AVGSES+VARS+SCOV(I)*.2D1 ENDIF ENDIF ENDIF 150 CONTINUE 140 FORMAT(1X,I4,F7.4,5X,F6.4,11X,F6.4,' - ',F6.4,15X, , 5X,F8.4,2(5X,F6.4),F17.4,' -',F7.4) 141 FORMAT(1X,I4,F7.4,5X,A6,11X,A6,' - ',A6,15X, , 5X,F8.4,F11.4,5X,A6,11X,A6,' - ',A6) 142 FORMAT(1X,I4,A7,5X,A6,11X,A6,' - ',A6,15X, , 5X,F8.4,2(5X,A6),11X,A6,' - ',A6) IF(AVGSPN.GT.0)THEN AVGT=AVGT/AVGSPN AVGSP=AVGSP/AVGSPN IF(AV1)THEN AVSESP=DSQRT(AVSESP)/AVSSPN AVST=DSQRT(AVST)/AVSSPN CI1=AVGSP-.196D1*AVSESP CI2=AVGSP+.196D1*AVSESP WRITE(7,160)AVGSP,AVSESP,CI1,CI2,AVGT,AVST, , AVGT-.196D1*AVST,AVGT+.196D1*AVST ELSE WRITE(7,161)AVGSP,(UNDEF,J=1,3),AVGT,(UNDEF,J=1,3) ENDIF ELSE WRITE(7,162)(UNDEF,J=1,8) ENDIF 160 FORMAT(/1X,130('-')/' MEAN',F7.4,5X,F6.4,F17.4,' -',F7.4, , 28X,2F11.4,F17.4,' -',F7.4/1X,130('-')) 161 FORMAT(/1X,130('-')/' MEAN',F7.4,A11,A17,' - ',A6, , 28X,F11.4,5X,A6,11X,A6,' - ',A6/1X,130('-')) 162 FORMAT(/1X,130('-')/' MEAN',A7,A11,A17,' - ',A6, , 28X,2(5X,A6),11X,A6,' - ',A6/1X,130('-')) C WRITE(7,170) 170 FORMAT(/8X,'==== Survival rate estimates between sampling ', ,'periods ==== ', ,'Interval Survival rate estimates per unit time (PHI**(1/T(i))'/ /' Period PHI SE(PHI) 95% Conf. ', ,'Interval COV(PHI(i,i-1)) Length PHI SE(PHI)', ,' 95% Conf. interval COV(PHI(i,i-1))'/1X,132('-')) AVGT=.0D0 AVST=.0D0 DO 180 I=2,NTIME-2 IF(CA(I+1,1)*NM(I)*(NM(I+1)+Z(I+1)).GT..0D0)THEN PHIT=.0D0 IF(S(I).GT..0D0)PHIT=S(I)**(.1D1/T(I)) AVGT=AVGT+PHIT IF(CA(I,1)*NM(I+1).GT..0D0)THEN CI1=S(I)-.196D1*SES(I) CI2=S(I)+.196D1*SES(I) TT=(.1D1-T(I))/T(I) XTMP=.0D0 IF(S(I-1).GT..0D0)XTMP=S(I-1)**((.1D1-T(I-1))/T(I-1))/T(I-1) TMP=.0D0 IF(S(I).GT..0D0)TMP=S(I)**TT/T(I) SEPHIT=SES(I)*TMP AVST=AVST+SEPHIT*SEPHIT IF(I.GT.2)THEN COVT=SCOV(I)*TMP*XTMP AVST=AVST+COVT*.2D1 WRITE(7,185)YR(I),S(I),SES(I),CI1,CI2,SCOV(I),T(I),PHIT, , SEPHIT,PHIT-.196D1*SEPHIT,PHIT+.196D1*SEPHIT,COVT ELSE WRITE(7,184)YR(I),S(I),SES(I),CI1,CI2,T(I),PHIT, , SEPHIT,PHIT-.196D1*SEPHIT,PHIT+.196D1*SEPHIT ENDIF ELSE WRITE(7,186)YR(I),S(I),(UNDEF,J=1,4),T(I),PHIT,(UNDEF,J=1,4) AV1=.FALSE. ENDIF ELSE WRITE(7,187)YR(I),(UNDEF,J=1,5),T(I),(UNDEF,J=1,5) AV1=.FALSE. ENDIF IF(MOD(YR(I),5).EQ.0)WRITE(7,140) 180 CONTINUE 184 FORMAT(1X,I4,F7.4,F7.4,4X,F11.4,' -',F7.4,21X, , 4X,F9.4,F8.4,F9.4,F9.4,' -',F10.4,F18.10) 185 FORMAT(1X,I4,F7.4,F7.4,4X,F11.4,' -',F7.4,F21.10, , 4X,F9.4,F8.4,F9.4,F9.4,' -',F10.4,F18.10) 186 FORMAT(1X,I4,F7.4,A7,A15,' - ',A6,9X,A12, , 4X,F9.4,F8.4,2A9,' -',A10,6X,A12) 187 FORMAT(1X,I4,A7,A7,A15,' - ',A6,9X,A12, , 4X,F9.4,A8,2A9,' -',A10,6X,A12) 188 FORMAT(1X,132('-')/' MEAN ',F6.4,F7.4,F15.4,' -', ,F7.4,34X,F8.4,2F9.4,' -',F10.4/1X,132('-')) 189 FORMAT(1X,132('-')/' MEAN',F7.4,A7,9X,A6,' - ', ,A6,35X,F7.4,2A9,' -',A10/1X,132('-')) 191 FORMAT(1X,132('-')/' MEAN ',A6,A7,9X,A6,' - ', ,A6,35X,A7,2A9,' -',A10/1X,132('-')) IF(AVGSN.GT.0)THEN AVGS=AVGS/AVGSN AVGT=AVGT/AVGSN IF(AV1)THEN AVGSES=DSQRT(AVGSES)/(NTIME-2) AVST=DSQRT(AVST)/(NTIME-2) CI1=AVGS-.196D1*AVGSES CI2=AVGS+.196D1*AVGSES WRITE(7,188)AVGS,AVGSES,CI1,CI2,AVGT,AVST, , AVGT-.196D1*AVST,AVGT+.196D1*AVST ELSE WRITE(7,189)AVGS,(UNDEF,J=1,3),AVGT,(UNDEF,J=1,3) ENDIF ELSE WRITE(7,191)(UNDEF,J=1,8) ENDIF C WRITE(7,190) 190 FORMAT(//' Period p SE(p) ', ,7X,'95% Conf. interval'/1X,62('-')) AV1=.TRUE. DO 200 I=2,NTIME-1 IF(Z(I)+CA(I,1).GT..0D0)THEN P(I)=CA(I,1)/(Z(I)+CA(I,1)) AVGPN=AVGPN+.1D1 AVGP=AVGP+P(I) IF(CA(I,1).GT..0D0)THEN VARP=P(I)*P(I)*(.1D1/CA(I,1)-.1D1/(CA(I,1)+Z(I))) SEP(I)=DSQRT(VARP) AVGSEP=AVGSEP+VARP CI1=P(I)-.196D1*SEP(I) CI2=P(I)+.196D1*SEP(I) WRITE(7,140)YR(I),P(I),SEP(I),CI1,CI2 ELSE WRITE(7,141)YR(I),P(I),(UNDEF,J=1,3) AV1=.FALSE. ENDIF ELSE WRITE(7,142)YR(I),(UNDEF,J=1,4) AV1=.FALSE. ENDIF IF(MOD(YR(I),5).EQ.0)WRITE(7,140) 200 CONTINUE IF(AVGPN.GT..0D0)THEN AVGP=AVGP/AVGPN IF(AV1)THEN AVGSEP=DSQRT(AVGSEP)/(NTIME-2) WRITE(7,205)AVGP,AVGSEP,AVGP-.196D1*AVGSEP, , AVGP+.196D1*AVGSEP ELSE WRITE(7,206)AVGP,(UNDEF,J=1,3) ENDIF ELSE WRITE(7,207)(UNDEF,J=1,4) ENDIF 205 FORMAT(1X,62('-')/' MEAN',F7.4,5X,F6.4,11X,F6.4,' - ',F6.4/ / 1X,62('-')) 206 FORMAT(1X,62('-')/' MEAN',F7.4,A11,11X,A6,' - ',A6/ / 1X,62('-')) 207 FORMAT(1X,62('-')/' MEAN',A7,A11,11X,A6,' - ',A6/1X,62('-')) C CALL Goodness-of-fit TEST ROUTINE FOR THIS MODEL CALL GOFM2(Z,NM,CHI2,IDF2) RETURN END C$CONTROL SEGMENT 'MODA' SUBROUTINE MODELA(T,NM,NN,NS,Z,ROW,AVGT1,AVGP,CHIA,IDFA,MODA,PRT, , CHIA2,IDFA2) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER YR,I,J,K,L,M,NTIME,NM1,NM2,JP1,IDF,IDFA,MODA,NP,NB,NPHI, , PRT,ITMP,IDFA2,IDIM PARAMETER (IDIM=50) LOGICAL VOK,VNFLG,AV1,AV2,AV3,PFLG,BFLG CHARACTER*133 LINE,UNDEF*6 DIMENSION NM(IDIM),NN(IDIM),NS(IDIM),Z(IDIM),ROW(IDIM),PHI(IDIM) DIMENSION COVN(IDIM,IDIM),AEST(IDIM),XN(IDIM),XM(IDIM) DIMENSION T(IDIM),FARR(1:IDIM,0:IDIM-2) COMMON /BLK0/NTIME,YR(IDIM),TITLE(16) EQUIVALENCE (FARR,COVN) DATA UNDEF/' Div/0'/ DATA AVG,AVST,SEFAVG,AVSDF2,SENAVG,AVSDN2,SEBAVG/7*.0D0/ DATA AVGN,AVGM,SEAVGM,SEAVGP,AVGB/5*.0D0/ DATA NP,NB,NPHI/3*0/,XM/IDIM*.0D0/,CNT9/.0D0/ DATA VOK,AV1,AV2/3*.TRUE./ C STATEMENT FUNCTION BEE(I)=FARR(I+1,I) C AVGT=.0D0 AVGT1=.0D0 AVGP=.0D0 WRITE(7,20)TITLE 20 FORMAT(/1X,16A8/) WRITE(7,10) 10 FORMAT(/' -- Model A -- Jolly-Seber capture-recapture model', ,' with both death and immigration'//' Div/0 Denotes estimates', /' which cannot be computed due to poor data (EG R(I)=0). Means', ,' do not include these time periods.') C C*****COMPUTE NUMBER ANIMALS MARKED IN SAMPLE I / NUMBER ANIMALS C CAUGHT IN SAMPLE I, AEST(I) C DO 30 I=2,NTIME AEST(I)=(NM(I)+.1D1)/(NN(I)+.1D1) 30 CONTINUE C C*****COMPUTE ESTIMATES OF SURV. RATES, PHI(I), AND POP. SIZE XN(I) C*****COMPUTE ESTIMATE OF TOTAL MARKED ANIMALS IN POPULATION AT C TIME I, XM(I) C WRITE(7,40) 40 FORMAT(/8X,'==== Survival rate estimates between sampling ', ,'periods ==== Interval', ,4X,'Survival rate estimates per unit time (PHI**(1/T(i))'/ /' Period PHI SE(PHI) ',8HSE'(PHI),' 95% Conf. ', ,'interval COV(PHI(i,i-1)) Length PHI SE(PHI)', ,' 95% Conf. interval COV(PHI(i,i-1))'/1X,132('-')) NM1=NTIME-1 NM2=NTIME-2 XN(1)=0 XM(1)=0 DO 50 I=2,NM1 VOK=.TRUE. XM(I)=(NS(I)+.1D1)*Z(I)/(ROW(I)+.1D1)+NM(I) XN(I)=DMAX1(XM(I)/AEST(I),NN(I)) IF(NS(I-1)*ROW(I-1).LT..1D1)THEN FARR(I,I-1)=-.5D55 IF(I.NE.2)WRITE(7,45)YR(I-1),(UNDEF,J=1,6),T(I-1), , (UNDEF,J=1,5) IF(I.EQ.2)WRITE(7,45)YR(I-1),(UNDEF,J=1,5),' ',T(I-1), , (UNDEF,J=1,4) ELSE PHI(I-1)=XM(I)/(NS(I-1)*Z(I-1)/ROW(I-1)+NS(I-1)) NPHI=NPHI+1 FARR(I,I-1)=XN(I)-PHI(I-1)*(XN(I-1)-NN(I-1)+NS(I-1)) PHIT=PHI(I-1)**(.1D1/T(I-1)) AVGT=AVGT+PHIT AVGT1=AVGT1+DMIN1(PHIT,.1D1) AVG=AVG+PHI(I-1) IF((Z(I)+NM(I))*ROW(I)*NS(I)*(Z(I-1)+NS(I-1)).GT.0)THEN C=.1D1/ROW(I-1) - .1D1/NS(I-1) C***** CALCULATE VARIANCE OF PHI VPHI2=XM(I)-NM(I) VPHI2=VPHI2*(VPHI2+NS(I))/(XM(I)*XM(I)) VPHI2=VPHI2*(.1D1/ROW(I) - .1D1/NS(I)) VPHI2=VPHI2+(XM(I-1)-NM(I-1))/ / (XM(I-1)-NM(I-1)+NS(I-1))*C VARPHI=PHI(I-1)*PHI(I-1)*(VPHI2+(.1D1-PHI(I-1))/XM(I)) VPHI2=VPHI2*PHI(I-1)*PHI(I-1) AVSDF2=AVSDF2+VPHI2 SEPHI2=DSQRT(VPHI2) SEFAVG=SEFAVG+VARPHI SEPHI=DSQRT(VARPHI) CI1=PHI(I-1)-.196D1*SEPHI CI2=PHI(I-1)+.196D1*SEPHI TT=(.1D1-T(I-1))/T(I-1) IF(PHI(I-1).LE..0D0)PHIT=.0D0 SEPHIT=SEPHI*(PHI(I-1)**TT/T(I-1)) AVST=AVST+SEPHIT*SEPHIT ELSE AV2=.FALSE. VOK=.FALSE. ENDIF IF(I.EQ.2)THEN IF(VOK)WRITE(7,60)YR(I-1),PHI(I-1),SEPHI,SEPHI2,CI1,CI2, , T(I-1),PHIT,SEPHIT,PHIT-.196D1*SEPHIT,PHIT+.196D1*SEPHIT IF(.NOT.VOK)WRITE(7,64)YR(I-1),PHI(I-1),(UNDEF,J=1,4),' ', , T(I-1),PHIT,(UNDEF,J=1,3) ELSE IF(NS(I-2)*ROW(I-2)*(Z(I-1)+NM(I-1)).GT..0D0)THEN COVPHI=-((PHI(I-2)*PHI(I-1)*(XM(I-1)-NM(I-1)))/ / XM(I-1))*C COVPHT=COVPHI*PHI(I-2)**((.1D1-T(I-2))/T(I-2))* * PHI(I-1)**TT /T(I-2)/T(I-1) AVST=AVST+COVPHT*.2D1 SEFAVG=SEFAVG+COVPHI*.2D1 AVSDF2=AVSDF2+COVPHI*.2D1 IF(VOK)WRITE(7,61)YR(I-1),PHI(I-1),SEPHI,SEPHI2,CI1,CI2, , COVPHI,T(I-1),PHIT,SEPHIT,PHIT-.196D1*SEPHIT, , PHIT+.196D1*SEPHIT,COVPHT IF(.NOT.VOK)WRITE(7,62)YR(I-1),PHI(I-1),(UNDEF,J=1,4), , COVPHI,T(I-1),PHIT,(UNDEF,J=1,3),COVPHT ELSE IF(VOK)WRITE(7,63)YR(I-1),PHI(I-1),SEPHI,SEPHI2,CI1,CI2, , UNDEF,T(I-1),PHIT,SEPHIT,PHIT-.196D1*SEPHIT, , PHIT+.196D1*SEPHIT,UNDEF IF(.NOT.VOK)WRITE(7,64)YR(I-1),PHI(I-1),(UNDEF,J=1,5), , T(I-1),PHIT,(UNDEF,J=1,4) ENDIF ENDIF ENDIF IF(MOD(YR(I-1),5).EQ.0)WRITE(7,60) 50 CONTINUE IF(NPHI.GT.0)THEN AVGT=AVGT/(NPHI) AVGT1=AVGT1/(NPHI) AVG=AVG/(NPHI) IF(AV2)THEN SEFAVG=DSQRT(SEFAVG)/(NM2) AVSDF2=DSQRT(AVSDF2)/(NM2) AVST=DSQRT(AVST)/(NM2) WRITE(7,70)AVG,SEFAVG,AVSDF2,AVG-.196D1*SEFAVG, ,AVG+.196D1*SEFAVG,AVGT,AVST,AVGT-.196D1*AVST,AVGT+.196D1*AVST ELSE WRITE(7,71)AVG,(UNDEF,J=1,4),AVGT,(UNDEF,J=1,3) ENDIF ELSE WRITE(7,72)(UNDEF,J=1,9) ENDIF 45 FORMAT(1X,I4,1X,A6,1X,A6,2X,A6,5X,A6,' - ',A6,11X,A6, , F11.4,2X,A6,2X,A6,4X,A6,' - ',A7,12X,A6) 60 FORMAT(1X,I4,F7.4,1X,F6.4,2X,F6.4,4X,F7.4,' -',F7.4,5X,12X, , 3X,F8.4,2F8.4,1X,F9.4,' -',F8.4,4X,F14.10) 61 FORMAT(1X,I4,F7.4,1X,F6.4,2X,F6.4,4X,F7.4,' -',F7.4,5X,F12.10, , 3X,F8.4,2F8.4,1X,F9.4,' -',F8.4,4X,F14.10) 62 FORMAT(1X,I4,F7.4,1X,A6,2X,A6,5X,A6,' - ',A6,5X,F12.10, , 3X,F8.4,F8.4,2X,A6,4X,A6,' - ',A7,4X,F14.10) 63 FORMAT(1X,I4,F7.4,1X,F6.4,2X,F6.4,4X,F7.4,' -',F7.4,5X,A12, , 3X,F8.4,2F8.4,1X,F9.4,' -',F8.4,4X,A14) 64 FORMAT(1X,I4,F7.4,1X,A6,2X,A6,5X,A6,' - ',A6,5X,A12, , 3X,F8.4,F8.4,2X,A6,4X,A6,' - ',A7,4X,A14) 70 FORMAT(1X,132('-')/' MEAN ',2(F6.4,1X),1X,F6.4,4X,F7.4,' - ', ,F6.4,28X,2F8.4,1X,F9.4,' -',F8.4/1X,132('-')) 71 FORMAT(1X,132('-')/' MEAN ',F6.4,1X,A6,1X,1X,A6,5X,A6,' - ', ,A6,28X,F8.4,2X,A6,4X,A6,' - ',A7/1X,132('-')) 72 FORMAT(1X,132('-')/' MEAN ',A6,1X,A6,1X,1X,A6,5X,A6,' - ', ,A6,28X,2X,A6,2X,A6,4X,A6,' - ',A7/1X,132('-')) C C*****COMPUTE EXPECTED NUMBER IN POPULATION AT TIME I WHICH FIRST C**** JOINED POPULATION BETWEEN TIMES J AND J+1,(FARR(I,J)) C XN(1)=XN(2) FARR(1,0)=XN(2) IF(NS(1)*ROW(1).GT.0.)THEN FARR(2,1)=XN(2)-PHI(1)*(XN(2)-NN(1)+NS(1)) IF(BEE(1).EQ..0D0)FARR(2,1)=XN(2)-PHI(1)*(XN(2)+1-NN(1)+NS(1)) ENDIF DO 160 I=0,NM2 DO 155 J=I+2,NM1 FARR(J,I)=-.5D55 IF(BEE(J-1).GE.-.5D54.AND.FARR(J-1,I).GE.-.5D54) ) FARR(J,I)=((XN(J)-BEE(J-1))/XN(J-1))*FARR(J-1,I) 155 CONTINUE 160 CONTINUE AV2=.TRUE. DO 200 I=2,NM1 DO 190 J=1,I-1 COVN(J,I)=.0D0 DO 180 K=0,J-1 IF(BEE(K).EQ..0D0.OR.FARR(I,K).LT.-.5D54.OR. * FARR(J,K).LT.-.5D54.OR.COVN(J,I).LT.-.5D54)THEN AV2=.FALSE. COVN(J,I)=-.5D55 GO TO 190 ELSE COVN(J,I)=COVN(J,I)+FARR(I,K)-FARR(J,K)*FARR(I,K)/BEE(K) ENDIF 180 CONTINUE IF(AV2)SENAVG=SENAVG+.2D1*COVN(J,I) 190 CONTINUE 200 CONTINUE C WRITE(7,120) 120 FORMAT(//' Period M',11X,'SE',1H','(M)', ,10X,'95% Conf. interval',13X,'N',11X,'SE(N) ',6HSE'(N),11X, ,'95% Conf. interval'/1X,132('-')) C AV1=.TRUE. AV3=.TRUE. DO 250 I=2,NM1 AVGM=AVGM+XM(I) AVGN=AVGN+XN(I) VNFLG=.TRUE. IF(ROW(I)*NS(I).GT..0D0)THEN VARM=(XM(I)-NM(I))*(XM(I)-NM(I)+NS(I))* * (.1D1/ROW(I)-.1D1/NS(I)) SEM=DSQRT(VARM) SEAVGM=SEAVGM+VARM IF((Z(I)*NM(I))*NM(I).GT..0D0)THEN VARN=XN(I)*(XN(I)-NN(I))*((XM(I)-NM(I)+NS(I))/XM(I)* * (.1D1/ROW(I)-.1D1/NS(I))+(.1D1-AEST(I))/NM(I)) SEN2=DSQRT(VARN) AVSDN2=AVSDN2+VARN SUM=.0D0 DO 125 J=0,I-1 IF(BEE(J).EQ..0D0.OR.FARR(I,J).LT.-.5D54)THEN VNFLG=.FALSE. ELSE SUM=SUM+FARR(I,J)*FARR(I,J)/BEE(J) ENDIF 125 CONTINUE IF(XN(I).GT.SUM)VARN=VARN+XN(I)-SUM SENAVG=SENAVG+VARN SEN=DSQRT(VARN) IF(VNFLG)THEN WRITE(7,130)YR(I),XM(I),SEM,XM(I)-.196D1*SEM,XM(I)+ + .196D1*SEM,XN(I),SEN,SEN2,XN(I)-.196D1*SEN,XN(I)+.196D1*SEN ELSE AV2=.FALSE. WRITE(7,131)YR(I),XM(I),SEM,XM(I)-.196D1*SEM,XM(I)+ + .196D1*SEM,XN(I),UNDEF,SEN2,UNDEF,UNDEF ENDIF ELSE AV2=.FALSE. AV3=.FALSE. WRITE(7,132)YR(I),XM(I),SEM,XM(I)-.196D1*SEM,XM(I)+ + .196D1*SEM,XN(I),(UNDEF,J=1,4) ENDIF ELSE AV1=.FALSE. AV2=.FALSE. WRITE(7,133)YR(I),XM(I),(UNDEF,J=1,3),XN(I),(UNDEF,J=1,4) ENDIF IF(MOD(YR(I),5).EQ.0)WRITE(7,130) 250 CONTINUE 130 FORMAT(1X,I4,F12.2,1X,F12.2,1X,F13.2,' - ',F11.2,9X, , F9.2,1X,2F12.2,1X,F12.2,' - ',F14.2) 131 FORMAT(1X,I4,F12.2,1X,F12.2,1X,F13.2,' - ',F11.2,9X, , F9.2,7X,A6,F12.2,7X,A6,' - ',8X,A6) 132 FORMAT(1X,I4,F12.2,1X,F12.2,1X,F13.2,' - ',F11.2,9X, , F9.2,1X,2(6X,A6),7X,A6,' - ',8X,A6) 133 FORMAT(1X,I4,F12.2,7X,A6,8X,A6,' - ',5X,A6,9X, , F9.2,1X,2(6X,A6),7X,A6,' - ',8X,A6) AVGM=AVGM/(NM2) AVGN=AVGN/(NM2) IF(AV1)THEN SEAVGM=DSQRT(SEAVGM)/(NM2) IF(AV3)THEN AVSDN2=DSQRT(AVSDN2)/(NM2) IF(AV2)THEN SENAVG=DSQRT(SENAVG)/(NM2) WRITE(7,260)AVGM,SEAVGM,AVGM-.196D1*SEAVGM,AVGM+.196D1* * SEAVGM,AVGN,SENAVG,AVSDN2,AVGN-.196D1*SENAVG,AVGN+.196D1*SENAVG ELSE WRITE(7,261)AVGM,SEAVGM,AVGM-.196D1*SEAVGM,AVGM+.196D1* * SEAVGM,AVGN,UNDEF,AVSDN2,UNDEF,UNDEF ENDIF ELSE WRITE(7,262)AVGM,SEAVGM,AVGM-.196D1*SEAVGM,AVGM+.196D1* * SEAVGM,AVGN,(UNDEF,J=1,4) ENDIF ELSE WRITE(7,263)AVGM,(UNDEF,J=1,3),AVGN,(UNDEF,J=1,4) ENDIF 260 FORMAT(1X,132('-')/' MEAN',F12.2,1X,F12.2,1X,F13.2,' - ', ,F11.2,9X,F9.2,1X,2F12.2,1X,F12.2,' - ',F14.2/1X,132('-')) 261 FORMAT(1X,132('-')/' MEAN',F12.2,1X,F12.2,1X,F13.2,' - ', ,F11.2,9X,F9.2,7X,A6,F12.2,7X,A6,' - ',8X,A6/1X,132('-')) 262 FORMAT(1X,132('-')/' MEAN',F12.2,1X,F12.2,1X,F13.2,' - ', ,F11.2,9X,F9.2,7X,A6,6X,A6,7X,A6,' - ',8X,A6/1X,132('-')) 263 FORMAT(1X,132('-')/' MEAN',F12.2,7X,A6,8X,A6,' - ', ,5X,A6,9X,F9.2,7X,A6,6X,A6,7X,A6,' - ',8X,A6/1X,132('-')) WRITE(7,80) 80 FORMAT(//' Period p SE(p)',5X, ,'95% Conf. interval',9X,'B',10X,'SE(B)',10X, ,'95% Conf. interval',7X,'COV(B(i),B(i-1))'/1X,132('-')) AV1=.TRUE. AV2=.TRUE. DO 300 I=2,NM1 WRITE(LINE,'(I5)')YR(I) IF(XM(I).GT..0D0)THEN P=NM(I)/XM(I) NP=NP+1 WRITE(LINE(6:),'(F11.4)')P AVGP=AVGP+P IF(ROW(I)*NS(I)*NM(I)*Z(I).GT..0D0)THEN VARP=P*P*(.1D1-P)*(.1D1-P)* * (.1D1/ROW(I)-.1D1/NS(I)+.1D1/NM(I)+.1D1/Z(I)) SEP=DSQRT(VARP) WRITE(LINE(17:),'(F11.4)')SEP WRITE(LINE(28:),'(F12.4)')P-.196D1*SEP WRITE(LINE(40:),'(3H - ,F6.4)')P+.196D1*SEP SEAVGP=SEAVGP+VARP ELSE WRITE(LINE(17:),'(A11)')UNDEF WRITE(LINE(28:),'(A12)')UNDEF WRITE(LINE(40:),'(3H - ,A6)')UNDEF AV1=.FALSE. ENDIF ELSE WRITE(LINE(6:),'(A11)')UNDEF WRITE(LINE(17:),'(A11)')UNDEF WRITE(LINE(28:),'(A12)')UNDEF WRITE(LINE(40:),'(3H - ,A6)')UNDEF AV1=.FALSE. ENDIF IF(I.EQ.NM1)GO TO 280 IF(NS(I)*ROW(I).GT..0D0)THEN B=BEE(I) AVGB=AVGB+B NB=NB+1 WRITE(LINE(49:),'(F14.2)')B C***** CALCULATE VARIANCE OF B IF(ROW(I-1)*NS(I-1)*NM(I)*NM(I+1)*ROW(I+1)*NS(I+1)* * (XM(I)-NM(I)+NS(I)).GT..0D0)THEN C=.1D1/ROW(I-1) - .1D1/NS(I-1) E=XM(I+1)-NM(I+1) F=(E+NS(I+1))/(XM(I+1)*XM(I+1)) G=(.1D1/ROW(I+1))-(.1D1/NS(I+1)) H=(XM(I)-NM(I))/(XM(I)-NM(I)+NS(I)) PP=((PHI(I)*NS(I)*(.1D1-AEST(I)))/AEST(I)) Q=(XN(I)-NN(I))*(XN(I+1)-B)*(.1D1-AEST(I))*(.1D1-PHI(I)) Q=Q/(XM(I)-NM(I)+NS(I)) R=XN(I+1)*(XN(I+1)-NN(I+1))*(.1D1-AEST(I+1))/NM(I+1) S=(PHI(I)*PHI(I))*XN(I)*(XN(I)-NN(I))*(.1D1-AEST(I))/NM(I) IF(Q.LT..0D0)Q=.0D0 IF(R.LT..0D0)R=.0D0 IF(S.LT..0D0)S=.0D0 VARBEE=(B*B*E*F*G)+(H*PP*PP*C)+Q+R+S SEB=DSQRT(VARBEE) WRITE(LINE(63:),'(F13.2)')SEB WRITE(LINE(76:),'(F14.2)')B-.196D1*SEB WRITE(LINE(90:),'(3H - ,F11.2)')B+.196D1*SEB SEBAVG=SEBAVG+VARBEE IF(I.GT.2.AND.I.LT.NM1)THEN COVB=-PHI(I)*(XN(I)-NN(I))*(.1D1-AEST(I))* 1 (BEE(I-1)*NS(I)/XM(I)*(1./ROW(I)-1./NS(I))+XN(I)/NM(I)) WRITE(LINE(104:),'(8X,F14.2)')COVB SEBAVG=SEBAVG+.2D1*COVB ENDIF ELSE WRITE(LINE(63:),'(A13)')UNDEF WRITE(LINE(76:),'(A14)')UNDEF WRITE(LINE(90:),'(3H - ,A11)')UNDEF IF(I.LT.NM2)WRITE(LINE(104:),'(8X,A14)')UNDEF AV2=.FALSE. ENDIF ELSE WRITE(LINE(49:),'(A14)')UNDEF WRITE(LINE(63:),'(A13)')UNDEF WRITE(LINE(76:),'(A14)')UNDEF WRITE(LINE(90:),'(3H - ,A11)')UNDEF IF(I.LT.NM2)WRITE(LINE(104:),'(8X,A14)')UNDEF AV2=.FALSE. ENDIF 280 WRITE(7,FMT='(A133)')LINE 300 IF(MOD(YR(I),5).EQ.0)WRITE(7,130) WRITE(7,FMT='(1X,132(1H-))') WRITE(LINE,'(1X,A4)')'MEAN' IF(NP.GT.0)THEN AVGP=AVGP/(NP) WRITE(LINE(6:),'(F11.4)')AVGP ELSE WRITE(LINE(6:),'(A11)')UNDEF ENDIF IF(AV1.AND.NP.GT.1)THEN SEAVGP=DSQRT(SEAVGP)/(NP-1) WRITE(LINE(17:),'(F11.4)')SEAVGP WRITE(LINE(28:),'(F12.4)')P-.196D1*SEAVGP WRITE(LINE(40:),'(3H - ,F6.4)')P+.196D1*SEAVGP ELSE WRITE(LINE(17:),'(A11)')UNDEF WRITE(LINE(28:),'(A12)')UNDEF WRITE(LINE(40:),'(3H - ,A6)')UNDEF ENDIF IF(NB.GT.0)THEN AVGB=AVGB/(NB) WRITE(LINE(49:),'(F14.2)')AVGB ELSE WRITE(LINE(49:),'(A14)')UNDEF ENDIF IF(AV2.AND.NB.GT.1)THEN SEBAVG=DSQRT(SEBAVG)/(NB-1) WRITE(LINE(63:),'(F13.2)')SEBAVG WRITE(LINE(76:),'(F14.2)')AVGB-.196D1*SEBAVG WRITE(LINE(90:),'(3H - ,F11.2)')AVGB+.196D1*SEBAVG ELSE WRITE(LINE(63:),'(A13)')UNDEF WRITE(LINE(76:),'(A14)')UNDEF WRITE(LINE(90:),'(3H - ,A11)')UNDEF ENDIF WRITE(7,FMT='(A133/1X,132(1H-))')LINE C C *** WRITE Covariance MATRIX IF(PRT.EQ.1)THEN WRITE(7,20)TITLE WRITE(7,10) WRITE(7,310) 310 FORMAT(/' Covariance(N(i),N(j)) (i=2,s-2 ; j=i+1,s-1)'/) K=6 DO 320 I=2,NM2 WRITE(7,*) WRITE(LINE,'(I5)')YR(I) DO 315 J=I+1,NM1 IF(COVN(I,J).LT.-.5D54)THEN WRITE(LINE(K:),'(A12)')UNDEF ELSE WRITE(LINE(K:),'(F12.2)')COVN(I,J) ENDIF K=K+12 IF(J.EQ.NM1.OR.K+12.GT.132)THEN K=6 WRITE(7,FMT='(A133)')LINE ENDIF 315 CONTINUE 320 CONTINUE ENDIF 340 WRITE(7,20)TITLE CALL GOF2(0,Z,Z,CHIA2,IDFA2,IDFA,CHIA) IF(IDFA.GT.0.AND.MODA.EQ.1)GO TO 360 WRITE(7,20)TITLE CALL GFIT(Z,ROW,NM,0,TMP,I) IF(IDFA.LE.0)THEN IDFA=I CHIA=TMP ENDIF 360 RETURN END C$CONTROL SEGMENT 'MODB' SUBROUTINE MODELB(T,M,N,NR,R,Z,STRT,EPS,MAXITR,PHIM1,CC,PRT) C C MODEL B: CONSTANT SURVIVAL RATE PER UNIT TIME, C TIME SPECIFIC CAPTURE PROBS. C C S = # OF SAMPLING PERIODS C FI = CONSTANT SURVIVAL RATE C T(I) = TIME UNITS BETWEEN SAMPLES I & I+1 C PHIM1(I) = SURVIVAL RATES FROM MODEL A (USED FOR C STARTING VALUE OF FI. C PEST(I) = CAPT. PROBS. FROM MODEL A (USED FOR C STARTING VALUE OF P(I). C 1-X(I) = Prob. THAT AN ANIMAL ALIVE JUST after C SAMPLE I IS SUBSEQUENTLY RECAPTURED. C STRT(I) = USER INPUT STARTING VALUES FOR FI & P(I) C NR(I) = # RELEASED FROM SAMPLE I C M(I) = # MARKED CAUGHT IN SAMPLE I C N(I) = # MARKED + UNMARKED CAUGHT IN SAMPLE I C R(I) = # CAUGHT IN SAMPLE I WHICH WERE LATER RECAPTURED C Z(I) = # CAUGHT before & after, BUT NOT IN SAMPLE I C G(I) = VECTOR OF FIRST ORDER PARTIAL DERATIVES C H(I,J) = MATRIX OF 2ND ORDER PARTIAL DERATIVES C V(I,J) = VARIANCE-Covariance MATRIX = INVERSE OF H C IMPLICIT DOUBLE PRECISION (A-Z) INTEGER PRT,ERR,YR,S,I,J,NITER,MAXITR,CC,TRY,IDIM PARAMETER (IDIM=50) DIMENSION T(IDIM),G(IDIM),M(IDIM),N(IDIM),NR(IDIM),R(IDIM), , Z(IDIM),STRT(IDIM),THETA(IDIM),THETA1(IDIM),H(IDIM,IDIM), , VCFIP(IDIM) COMMON /BLK0/S,YR(IDIM),TITLE(16) COMMON /BLK4/X(IDIM),P(IDIM),XD(IDIM),PD(1) DATA TRY/1/ C C STATEMENT FUNCTIONS C Q(I)=.1D1-P(I) PHI(I)=DSIGN(DABS(FI)**ABS(T(I)),FI) CIL(I)=THETA(I)-.196D1*DSQRT(H(I,I)) CIU(I)=THETA(I)+.196D1*DSQRT(H(I,I)) CC=0 C DO 5 I=1,S DO 5 J=1,S 5 H(J,I)=.0D0 C OUTPUT HEADER INFO WRITE(7,20)CHAR(12),TITLE 20 FORMAT(A1,16A8/ //' Model B - Constant survival rate per unit time, ', 2'time-specific capture probabilities.'/ 3/' Parameters: PHI = Survival rate per unit time'/ 4 ' p(i) = Capture probability in period i'/) 30 NITER=0 C IF(STRT(1).EQ..0D0)TRY=3 C C COMPUTE STARTING VALUES FROM MODEL A. C IF(TRY.EQ.2.OR.TRY.EQ.3)THEN IF(TRY.EQ.3)STRT(1)=PHIM1 IF(STRT(1).LT..01.OR.STRT(1).GT..99)STRT(1)=.5 PSUM=.0D0 DO 50 I=2,S-1 IF(M(I)+Z(I).GT..0D0)THEN PEST=M(I)/(M(I)+(NR(I)+.1D1)/(R(I)+.1D1)*Z(I)) ELSE PEST=.01 ENDIF IF(PEST.LT..01)PEST=.01 IF(PEST.GT..99)PEST=.99 PSUM=PSUM+PEST STRT(I)=PEST 50 CONTINUE STRT(S)=PSUM/DBLE(FLOAT(S-2)) ENDIF C C SET STARTING VALUES FROM INPUT ARRAY (STRT) C FI=STRT(1) DO 70 I=2,S 70 P(I)=STRT(I) C C COMPUTE PROB OF NO-RECAPTURE - X(I) C 80 X(S)=.1D1 P(1)=P(S) DO 90 I=1,S-1 J=S-I X(J)=.1D1-PHI(J)*(.1D1-Q(J+1)*X(J+1)) 90 THETA(J+1)=X(J) THETA(1)=FI C C OUTPUT STARTING VALUE OF THETA C WRITE(7,100)FI,(P(I),I=2,S) 100 FORMAT(/' Starting value of PHI, p(i):'/(10(1X,F8.4))) C C COMPUTE PARAMETER VECTOR 'THETA' C 110 FI=THETA(1) X(1)=THETA(2) DO 120 I=2,S-1 X(I)=THETA(I+1) 120 P(I)=((.1D1-X(I-1))/PHI(I-1)-(.1D1-X(I)))/X(I) P(S)=(.1D1-X(S-1))/PHI(S-1) C C COMPUTE GRADIANT VECTOR 'G' AND C VARIANCE-Covariance MATRIX 'V' C SUM=.0D0 SUM2=.0D0 DO 130 I=2,S-1 SUM=SUM+T(I-1)/X(I)*(Z(I)/Q(I)-M(I)*(.1D1-X(I))/P(I)) SUM2=SUM2+T(I-1)**2/X(I)**2*(.1D1-Q(I)*X(I))* * (Z(I)/Q(I)**2+M(I)/P(I)**2*(.1D1-X(I))) 130 CONTINUE G(1)=SUM/FI H(1,1)=SUM2/FI**2 H(1,2)=T(1)/(FI*X(2)**2*PHI(1))* * (Z(2)/Q(2)**2+M(2)*(.1D1-X(2))/P(2)**2) DO 150 I=2,S G(I)=(Q(I-1)*M(I-1)/P(I-1)-Z(I-1)+NR(I-1)-R(I-1))/X(I-1)- - (Q(I)*M(I)/P(I)-Z(I))/(PHI(I-1)*Q(I)*X(I)) IF(I.GT.2)H(1,I)=-T(I-2)*M(I-1)*(.1D1-Q(I-1)*X(I-1))/ / (FI*P(I-1)**2*X(I-1)**2)+T(I-1)/(FI*X(I)**2*PHI(I-1))* * (Z(I)/Q(I)**2+M(I)*(.1D1-X(I))/P(I)**2) H(I,1)=H(1,I) H(I,I)=(M(I-1)/P(I-1)**2-M(I-1)-Z(I-1)+NR(I-1)-R(I-1))/ / X(I-1)**2+(M(I)/P(I)**2+Z(I)/Q(I)**2)/ / (PHI(I-1)*X(I))**2 IF(I.EQ.S)GO TO 150 H(I,I+1)=-M(I)/(PHI(I-1)*(P(I)*X(I))**2) H(I+1,I)=H(I,I+1) 150 CONTINUE C C COMPUTE NEXT THETA: 'THETA1' C IN MATRIX NOTATION... C THETA1 = THETA + V * G C 160 CALL MATINV(S,H) CALL MMULT(H,G,THETA1,S,S,1,IDIM,IDIM,1) CALL MATADD(THETA1,THETA,THETA1,S) C C CHECK FOR CONVERGENCE C CONVERGES WHEN ABS(THETA1(I)-THETA(I)) < EPS C FOR I=1,2,...S C ERR=0 DO 170 I=1,S IF(H(I,I).LT..0D0)ERR=1 IF(DABS(THETA(I)-THETA1(I)).GT.EPS)GO TO 180 170 CONTINUE IF(ERR.EQ.0)GO TO 240 IF(ERR.NE.0)GO TO 210 C C NOT CONVERGED... INCREASE # ITERATIONS AND C COMPUTE NEW THETA C 180 NITER=NITER+1 IF(NITER.GT.MAXITR)GO TO 210 IF(THETA1(1).LT..0D0)GO TO 210 DO 200 I=1,S IF(DABS(THETA1(I)).GT..1D2)GO TO 210 200 THETA(I)=THETA1(I) GO TO 110 C 210 WRITE(7,220)TRY,NITER,(THETA(I),I=1,S) 220 FORMAT(/' (',I2,')Convergence failure after ',I6,' iterations'/ / ' THETA = ',(1X,10G12.6)) WRITE(7,230)(THETA1(I),I=1,S) 230 FORMAT(' THETA1 = ',(1X,10G12.6)) C----------------------Try different starting values ------------ C TRY=1:User input starting values (or model D if no input) C TRY=2:Model D starting value for PHI & Model A values for P(i) C TRY=3:Model A starting values for PHI & P(i) C TRY=4-7:Since there is no estimate of P(s) from model A, the mean C of the other P(i) was used in TRY=2&3. If the real value of C P(s) is not close to the mean, model B may not converge. C Hopefully, one of these values (.01,.31,.61,.91) will be C close enough to P(s) to cause convergence. C----------------------------------------------------------------- TRY=TRY+1 IF(TRY.LE.3)GO TO 30 IF(TRY.EQ.4)THEN STRT(S)=.01 ELSE STRT(S)=STRT(S)+.3 ENDIF IF(STRT(S).LT..1D1)GO TO 30 RETURN C C ESTIMATES CONVERGED... PRINT RESULTS C 240 WRITE(7,250)NITER 250 FORMAT(/' Final values after ',I6,' iterations:'/) CC=1 WRITE(7,260) 260 FORMAT(/30X,'Standard',12X,'95%'/' Parameter ', ,'Estimate Variance Error Confidence interval',6X,'Covar', ,'iance(p(i),PHI)',14X,'i T(i) PHI**T(i)'/1X,132('-')) WRITE(7,270)THETA(1),H(1,1),DSQRT(H(1,1)),CIL(1),CIU(1) 270 FORMAT(/' PHI',5X,F7.4,1X,F11.8,2X,F7.4,4X,F7.4,' - ',F7.4/) P(1)=.1D30 AVGP=.0D0 AVGVP=.0D0 DO 290 I=2,S VARP=H(I,I)/(PHI(I-1))**2+ + T(I-1)**2*(.1D1-Q(I)*X(I))**2/FI**2*H(1,1)+ + .2D1*T(I-1)*(.1D1-Q(I)*X(I))/(FI*PHI(I-1))*H(1,I) VCFIP(I)=-T(I-1)*(.1D1-Q(I)*X(I))/FI*H(1,1)- - H(1,I)/PHI(I-1) IF(I.EQ.S)GO TO 280 VARP=(VARP+Q(I)**2*H(I+1,I+1)- - .2D1*Q(I)/PHI(I-1)*H(I,I+1)- - .2D1*Q(I)*T(I-1)*(.1D1-Q(I)*X(I))/FI*H(1,I+1))/X(I)**2 VCFIP(I)=(VCFIP(I)+Q(I)*H(1,I+1))/X(I) 280 C1=P(I)-.196D1*DSQRT(VARP) C2=P(I)+.196D1*DSQRT(VARP) WRITE(7,310)YR(I),P(I),VARP,DSQRT(VARP),C1,C2,VCFIP(I), , I-1,T(I-1),PHI(I-1) AVGP=AVGP+P(I) AVGVP=AVGVP+VARP 290 IF(MOD(YR(I),5).EQ.0)WRITE(7,300) 300 FORMAT(' ') 310 FORMAT(' p(',I4,') ',F7.4,F12.8,2X,F7.4,4X,F7.4,' - ', ,F7.4,11X,F12.8,18X,I2,2X,F8.4,4X,F8.4) AVGP=AVGP/DBLE(FLOAT(S-1)) AVGVP=AVGVP/DBLE(FLOAT(S-1)) C1=AVGP-.196D1*DSQRT(AVGVP) C2=AVGP+.196D1*DSQRT(AVGVP) WRITE(7,320)AVGP,AVGVP,DSQRT(AVGVP),C1,C2 320 FORMAT(1X,58('-')/' p(MEAN) ',F7.4,F12.8,2X,F7.4,4X,F7.4,' - ', ,F7.4/1X,58('-')) C C CALL SUBROUTINE TO PRINT VAR-COV MATRICES CALL VARCVB(S,M,N,Z,FI,P,T,X,H,VCFIP,PRT) C RETURN END C$CONTROL SEGMENT 'VCB' SUBROUTINE MMULT(A,B,C,I1,I2,I3,N1,N2,N3) C C MATRIX MULTIPLICATION SUBROUTINE C A = N1 X N2 MATRIX C B = N2 X N3 MATRIX C C = A * B C DOUBLE PRECISION A(N1,N2),B(N2,N3),C(N1,N3) DO 10 I=1,I1 DO 10 J=1,I3 C(I,J)=.0D0 DO 10 K=1,I2 10 C(I,J)=C(I,J)+A(I,K)*B(K,J) RETURN END C$CONTROL SEGMENT 'VCB' SUBROUTINE MATADD(A,B,C,N) C C MATRIX ADDITION SUBROUTINE C A & B = IDIM X 1 MATRICES C C = A * B C PARAMETER (IDIM=50) DOUBLE PRECISION A(IDIM),B(IDIM),C(IDIM) DO 10 I=1,N 10 C(I)=A(I)+B(I) RETURN END C$CONTROL SEGMENT 'VCB' SUBROUTINE MATINV (M,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C---------------------------------------------------------- C INVERSE AND DETERMINANT OF A BY THE GAUSS-JORDAN METHOD. C M IS THE ORDER OF THE SQUARE MATRIX, A. C A-INVERSE REPLACES A. C DETERMINANT OF A IS PLACED IN DET. C COOLEY AND LOHNES (1971:63) C----------------------------------------------------------- PARAMETER (IDIM=50) DIMENSION A(IDIM,IDIM),IPVT(IDIM),PVT(IDIM),IND(IDIM,2) DET=.1D1 DO 10 J=1,M 10 IPVT(J)=0 DO 100 I=1,M C SEARCH FOR THE PIVOT ELEMENT AMAX=.0D0 DO 40 J=1,M IF(IPVT(J).EQ.1)GO TO 40 DO 30 K=1,M IF(IPVT(K)-1)20,30,130 20 IF(DABS(AMAX).GE.DABS(A(J,K)))GO TO 30 IROW=J ICOL=K AMAX=A(J,K) 30 CONTINUE 40 CONTINUE IPVT(ICOL)=IPVT(ICOL)+1 C INTERCHANGE THE ROWS TO PUT THE PIVOT ELEMENT ON THE DIAGONAL IF(IROW.EQ.ICOL)GO TO 60 DET=-DET DO 50 L=1,M SWAP=A(IROW,L) A(IROW,L)=A(ICOL,L) 50 A(ICOL,L)=SWAP 60 IND(I,1)=IROW IND(I,2)=ICOL PVT(I)=A(ICOL,ICOL) C DIVIDE THE PIVOT ROW BY THE PIVOT ELEMENT A(ICOL,ICOL)=.1D1 DO 70 L=1,M A(ICOL,L)=A(ICOL,L)/PVT(I) 70 CONTINUE C REDUCE THE NON-PIVOT ROWS DO 100 L1=1,M IF(L1.EQ.ICOL)GO TO 90 SWAP=A(L1,ICOL) A(L1,ICOL)=.0D0 IF(SWAP.LT..1D-30.AND.SWAP.GT.-.1D-30)SWAP=.0D0 DO 80 L=1,M A(L1,L)=A(L1,L)-A(ICOL,L)*SWAP 80 CONTINUE 90 CONTINUE 100 CONTINUE C INTERCHANGE THE COLUMNS DO 120 I=1,M L=M+1-I IF(IND(L,1).EQ.IND(L,2))GO TO 120 IROW=IND(L,1) ICOL=IND(L,2) DO 110 K=1,M SWAP=A(K,IROW) A(K,IROW)=A(K,ICOL) A(K,ICOL)=SWAP 110 CONTINUE 120 CONTINUE 130 CONTINUE RETURN END C$CONTROL SEGMENT 'MODD' SUBROUTINE MODELD(T,M,N,NR,R,Z,STRT,EPS,MAXITR,PHIM1,PEST,CC,PRT) C C MODEL D: CONSTANT SURVIVAL RATE PER UNIT TIME, C AND CONSTANT CAPTURE PROBS. C C S = # OF SMAPLING PERIODS C FI = CONSTANT SURVIVAL RATE C T(I) = TIME UNITS BETWEEN SAMPLES I & I+1 C 1-X(I) = Prob. THAT AN ANIMAL ALIVE JUST after C SAMPLE I IS SUBSEQUENTLY RECAPTURED. C NR(I) = # RELEASED FROM SAMPLE I C M(I) = # MARKED CAUGHT IN SAMPLE I C N(I) = # MARKED + UNMARKED CAUGHT IN SAMPLE I C R(I) = # CAUGHT IN SAMPLE I WHICH WERE LATER RECAPTURED C Z(I) = # CAUGHT before & after, BUT NOT IN SAMPLE I C G(I) = VECTOR OF FIRST ORDER PARTIAL DERATIVES C H(I,J) = MATRIX OF 2ND ORDER PARTIAL DERATIVES C V(I,J) = VARIANCE-Covariance MATRIX = INVERSE OF H C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) INTEGER S,CC,YR,PRT DOUBLE PRECISION T(IDIM),M(IDIM),N(IDIM),NR(IDIM) DIMENSION THETA(2),THETA1(2),G(2),V(2,2),W(IDIM),GMA(IDIM), , ALFA(IDIM,IDIM),R(IDIM),Z(IDIM),H(2,2),STRT(IDIM) EQUIVALENCE (V(1,1),H(1,1)) COMMON /BLK0/S,YR(IDIM),TITLE(16) COMMON /BLK4/XB(IDIM),PB(IDIM),X(IDIM),P C C STATEMENT FUNCTIONS C PHI(I)=FI**T(I) C C OUTPUT HEADER INFO C CC=0 WRITE(7,20)CHAR(12),TITLE 20 FORMAT(A1,16A8/ //' Model D - Constant survival rate per unit time, ', 2'constant capture probability.'/ 3/' Parameters: PHI = Estimate of survival rate'/ 4 ' p = Capture probability '/) WRITE(7,30) 30 FORMAT(/' Definitions:'/ 1' THETA = Vector of parameters'/ 2' 1-X(i) = Probability that an animal alive just after', 3' period i is subsequently recaptured'/ 4' T(i) = Time units between periods i and i+1'/) NITER=0 C C SET UP STARTING VALUES FOR FI,P C IF(STRT(1).GT..0D0)GO TO 60 C C COMPUTE STARTING VALUES FROM MODEL A. C T(S) = FLAG TO INDICATE UNEQUAL TIME PERIODS C FI=PHIM1 IF(FI.LT..01.OR.FI.GT..99)FI=.5 P=PEST IF(P.LT..01.OR.P.GT..99)P=.5 GO TO 70 C C SET STARTING VALUES FROM INPUT ARRAY (STRT) C 60 FI=STRT(1) P=STRT(2) 70 THETA(1)=FI THETA(2)=P C WRITE(7,80)FI,P 80 FORMAT(/' Starting values of THETA :',5X,2F12.4) C C COMPUTE X(I) C 90 Q=.1D1-P X(S)=.1D1 DO 110 I=S-1,1,-1 ALFA(I,I)=.1D1 IF(I.LE.1)GO TO 110 DO 100 K=I,2,-1 100 ALFA(K-1,I)=ALFA(K,I)*Q*PHI(K) 110 X(I)=.1D1-PHI(I)*(.1D1-Q*X(I+1)) C C COMPUTE GRADIANT VECTOR & 2ND DERATIVES C W(1)=.0D0 G(1)=.0D0 G(2)=.0D0 H(1,1)=.0D0 H(1,2)=.0D0 H(2,2)=.0D0 GMA(1)=.0D0 DO 160 I=2,S GMA(I)=PHI(I-1)*(Q*GMA(I-1)+(NR(I-1)-R(I-1))/X(I-1)) G(1)=G(1)+T(I-1)*(M(I)+Z(I)-GMA(I)*(.1D1-Q*X(I))) G(2)=G(2)+M(I)/P-Z(I)/Q-X(I)*GMA(I) W(I)=PHI(I-1)**2*(Q**2*W(I-1)+NR(I-1)/X(I-1)) H(1,1)=H(1,1)+ T(I-1)**2 *(.1D1-Q*X(I))*(GMA(I)+ + (.1D1-Q*X(I))*W(I)) SUM1=.0D0 DO 130 J=I,S 130 SUM1=SUM1+ALFA(I-1,J-1)*X(J) H(1,2)=H(1,2)+T(I-1)*((.1D1-Q*X(I))*W(I)+GMA(I))*SUM1 H(2,2)=H(2,2)+M(I)/P**2+Z(I)/Q**2+X(I)**2*W(I) IF(I.LE.2)GO TO 150 SUM2=.0D0 SUM3=.0D0 DO 140 J=2,I-1 SUM2=SUM2+PHI(J)*ALFA(J,I-1)*(Q*X(J)*W(J)-GMA(J)) 140 SUM3=SUM3+T(J-1)*ALFA(J-1,I-1)*(GMA(J)+(.1D1-Q*X(J))*W(J)) H(1,1)=H(1,1)+.2D1*T(I-1)*(.1D1-Q*X(I))*SUM3 H(1,2)=H(1,2)+T(I-1)*(.1D1-Q*X(I))*SUM2 H(2,2)=H(2,2)+.2D1*X(I)*SUM2 150 H(2,1)=H(1,2) 160 CONTINUE G(1)=G(1)/FI H(1,1)=H(1,1)/FI**2 H(1,2)=H(1,2)/FI C C COMPUTE NEW THETA = THETA + V * G C CALL MATI22(H) CALL MMULT(H,G,THETA1,2,2,1,2,2,1) CALL MATADD(THETA1,THETA,THETA1,2) C C CHECK FOR CONVERGENCE C IF(THETA1(1).LT..0D0.OR.THETA1(2).LT..0D0)GO TO 190 IF(THETA1(1).GT..1D2.OR.THETA1(2).GT..1D2)GO TO 190 IF(DABS(THETA(1)-THETA1(1)).GT.EPS.OR. . DABS(THETA(2)-THETA1(2)).GT.EPS)GO TO 180 GO TO 220 C C NOT CONVERGED... INC # ITERATIONS, COMPUTE NEW THETA C 180 NITER=NITER+1 IF(THETA1(1).LT..0D0)GO TO 190 IF(NITER.GT.MAXITR)GO TO 190 THETA(1)=THETA1(1) THETA(2)=THETA1(2) FI=THETA(1) P=THETA(2) GO TO 90 C 190 WRITE(7,200)NITER,THETA(1),THETA(2),THETA1(1),THETA1(2) 200 FORMAT(/' Convergence failure after ',I6,' iterations, ', , ' THETA = ',2(1X,G12.6),', THETA1 = ',2(1X,G12.6)) RETURN C C ESTIMATES CONVERGED... PRINT RESULTS C 220 WRITE(7,230)NITER 230 FORMAT(/' Final values after ',I6,' iterations:'/) CC=1 C1=FI-.196D1*DSQRT(V(1,1)) C2=FI+.196D1*DSQRT(V(1,1)) C3=P-.196D1*DSQRT(V(2,2)) C4=P+.196D1*DSQRT(V(2,2)) WRITE(7,240)FI,V(1,1),DSQRT(V(1,1)),C1,C2,V(1,2), , P,V(2,2),DSQRT(V(2,2)),C3,C4 240 FORMAT(/38X,'Standard'/' Parameter',17X,'Variance error',8X, ,'95% confidence interval Covariance(P,PHI)'/1X,97('-')/ //6X,'PHI',4X,F7.4,3X,F12.8,2X,F7.4,11X,F7.4,' - ',F7.4,6X,F15.8/ //6X,'p',6X,F7.4,3X,F12.8,2X,F7.4,11X,F7.4,' - ',F7.4/1X,75('-')) IF(T(S).NE.99.)GO TO 280 WRITE(7,250) 250 FORMAT(/' Period',9X,'T(i)',8X,'PHI**T(i)'/1X,36('-')) DO 260 I=1,S-1 WRITE(7,270)YR(I),T(I),PHI(I) 260 IF(MOD(YR(I),5).EQ.0)WRITE(7,270) 270 FORMAT(2X,I4,5X,F7.4,6X,F10.4) 280 IF(STRT(3).GT..0D0)GO TO 300 STRT(1)=FI DO 290 I=2,S 290 STRT(I)=P 300 CALL VARCVD(FI,P,M,N,NR,Z,T,X,V,PRT) RETURN END C$CONTROL SEGMENT 'MODD' SUBROUTINE MATI22(X) C C INVERT A 2 X 2 SYMETRIC MATRIX C DOUBLE PRECISION DET,A,X(2,2) DET=X(1,1)*X(2,2)-X(1,2)*X(2,1) A=X(2,2)/DET X(2,2)=X(1,1)/DET X(1,2)=-X(1,2)/DET X(2,1)=X(1,2) X(1,1)=A RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE TESTS(M,N,NS,R,Z,A1CC,M2CC,DCC,BCC,CHIA,IDFA, , CHIAP,IDFAP,CHI2,IDF2,CHIA2,IDFA2) C C COMPUTE Chi-square TEST OF MODEL D VS MODEL B C AND TEST OF MODEL B VS MODEL A C C VARIABLES: C DCC=MODEL D COMPLETION CODE (0-SUCC, 1-FAILURE) C BCC=MODEL B COMPLETION CODE C M(I),N(I),NS(I),R(I),Z(I)...DEFINED IN MODEL B C XB(I)=X ARRAY FROM MODEL B C XD(I)=X ARRAY FROM MODEL D C PB(I)=P ARRAY FROM MODEL B C PD =P FROM MODEL D C IMPLICIT DOUBLE PRECISION (A-Z) INTEGER YR,I,J,S,DCC,BCC,IDFA,IDFB,IDFD,DBL1DF,DBL2DF,DAL1DF, , DAL2DF,BAL1DF,BAL2DF,IDFAP,IDF2,A1CC,M2CC,IDFA2,IDIM PARAMETER (IDIM=50) DIMENSION N(IDIM),M(IDIM),NS(IDIM),R(IDIM),Z(IDIM),T(6,IDIM) COMMON /BLK0/S,YR(IDIM),TITLE(16) COMMON /BLK4/XB(IDIM),PB(IDIM),XD(IDIM),PD(1) DATA DAT1,DAT2,BAT1,BAT2,DAL1,DAL2,BAL1,BAL2,DBL1,DBL2, ,DBL1DF,DBL2DF,DAL1DF,DAL2DF,BAL1DF,BAL2DF/10*.0D0,6*0/ DATA T/IDIM*.0D0,IDIM*.0D0,IDIM*.0D0,IDIM*.0D0,IDIM*.0D0, , IDIM*.0D0/ C WRITE(7,10)CHAR(12),TITLE 10 FORMAT(A1,16A8) IF(DCC.LE.0.AND.BCC.LE.0)GO TO 95 IF(DCC.LE.0.OR.BCC.LE.0)GO TO 20 CALL TDBL1(S,R,NS,XD,XB,DBL1,DBL1DF,5,T) CALL TDBL2(S,M,Z,PD,PB,XD,XB,DBL2,DBL2DF,6,T) DBL1DF=S-2-DBL1DF-DBL2DF C DBL1=-.2D1*DBL1 C DBL2=-.2D1*DBL2 C CHI1=DBL1+DBL2 C IF(DCC.LE.0.OR.BCC.LE.0)GO TO 60 C WRITE(7,40) C WRITE(7,140)DBL1,DBL2,CHI1,DBL1DF,CHIPRB(DBL1DF,CHI1) C 40 FORMAT(//' TEST OF MODEL D VS. MODEL B:') C 20 IF(DCC.LE.0)GO TO 30 CALL TDA1(S,R,NS,XD,DAL1,DAT1,DAL1DF,1,T) CALL TDA2(S,M,Z,PD,XD,DAL2,DAT2,DAL2DF,2,T) C DAL1=-.2D1*DAL1 C DAL2=-.2D1*DAL2 DAL1DF=2*S-5 -DAL1DF-DAL2DF IF(DCC.LE.0)GO TO 30 C WRITE(7,80) C WRITE(7,140)DAL1,DAL2,DAL1+DAL2,DAL1DF, C , CHIPRB(DAL1DF,DAL1+DAL2) C 80 FORMAT(//' TEST OF MODEL D VS. MODEL A:') 30 IF(BCC.LE.0)GO TO 40 CALL TDA1(S,R,NS,XB,BAL1,BAT1,BAL1DF,3,T) CALL TDA2(S,M,Z,PB,XB,BAL2,BAT2,BAL2DF,4,T) C BAL1=-.2D1*BAL1 C BAL2=-.2D1*BAL2 BAL1DF=S-3 -BAL1DF-BAL2DF IF(BCC.LE.0)GO TO 40 C WRITE(7,120) C WRITE(7,140)BAL1,BAL2,BAL1+BAL2,BAL1DF, C , CHIPRB(BAL1DF,BAL1+BAL2) C 120 FORMAT(//' TEST OF MODEL B VS. MODEL A:') C 140 FORMAT( ' ============================'/ C / ' Component-1 Chi-square = ',F12.4/ C / ' Component-2 Chi-square = ',F12.4/ C / ' TOTAL Chi-square = ',F12.4/ C / ' Degrees of freedom = ',F8.0/ C / ' PROBABILITY = ',F12.4) 40 WRITE(7,50) 50 FORMAT(/14X,'Test of model D vs model A',8X, / 'Test of model B vs model A',8X,'Test of model D vs model B'/ //' Per ',2(' T1 T2 '), , ' T1 T2 '/ / ' ',2(' ------- -------'), , ' ------- -------') DO 60 I=1,S-1 60 WRITE(7,70)YR(I),(T(J,I),J=1,6) 70 FORMAT(1X,I5,6F17.4) WRITE(7,80)(T(J,IDIM),J=1,6) 80 FORMAT(/' Total',6F17.4) DAT1=T(1,IDIM)+T(2,IDIM) BAT1=T(3,IDIM)+T(4,IDIM) DBL1=T(5,IDIM)+T(6,IDIM) WRITE(7,90)DAT1,BAT1,DBL1,DAL1DF,BAL1DF,DBL1DF, , CHIPRB(DAL1DF,DAT1),CHIPRB(BAL1DF,BAT1),CHIPRB(DBL1DF,DBL1) 90 FORMAT(/3(5X,'Total chi-square = ',F12.4)/ / 3(5X,'Degrees of freedom = ',I8,4X)/ / 3(5X,'Probability = ',F12.4)) 95 CHIB=CHIA+BAT1 IDFB=BAL1DF+IDFA CHID=CHIA+DAT1 IDFD=DAL1DF+IDFA IF(A1CC.GT.0)CALL TSTA1(R,M,N,NS,S) IF(M2CC.GT.0)WRITE(7,97)CHIA2,IDFA2,CHIPRB(IDFA2,CHIA2) 97 FORMAT(/' Likelihood-Ratio test of model A versus model 2:'/ //' Hypothesis: Survival rate of newly captured animals = ', ,' survival rate of previously captured animals.'/ //15X,'Chi-square = ',F10.4,' with ',I3,' degrees of freedom ', ,'(Probability = ',F8.4,')') WRITE(7,100) WRITE(7,110)CHIA,IDFA,CHIPRB(IDFA,CHIA) IF(A1CC.GT.0)WRITE(7,120)CHIAP,IDFAP,CHIPRB(IDFAP,CHIAP) IF(M2CC.GT.0)WRITE(7,130)CHI2,IDF2,CHIPRB(IDF2,CHI2) IF(IDFA.LE.0)GO TO 160 IF(IDFB.GT.IDFA)WRITE(7,140)CHIB,IDFB,CHIPRB(IDFB,CHIB) IF(IDFD.GT.IDFA)WRITE(7,150)CHID,IDFD,CHIPRB(IDFD,CHID) 100 FORMAT(/' Goodness-of-fit tests:'// / ' Model Chi-square DF Probability'/ / ' ----- ---------- -- -----------') 110 FORMAT(' A ',F9.4,I8,4X,F8.4) 120 FORMAT(/9H A' ,F9.4,I8,4X,F8.4) 130 FORMAT(/' 2 ',F9.4,I8,4X,F8.4) 140 FORMAT(/' B ',F9.4,I8,4X,F8.4) 150 FORMAT(/' D ',F9.4,I8,4X,F8.4) 160 RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE TSTA1(R,M,N,NS,NSAMP) C C TEST OF MODEL A' VS MODEL A C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IDIM=50) DOUBLE PRECISION R(IDIM),M(IDIM),N(IDIM),NS(IDIM),T(IDIM),C(IDIM) REAL ZPROB,Z,Z1 WRITE(7,10) 10 FORMAT(/' Test of Model A versus A',1H',' : see Pollock (1974:', ,'Biometrics 30: 77-87)'//10X,'i',10X,'P1(i)',10X,'P2(i)',10X, ,'Z(i)',10X,'Pr(Z(i))'/5X,68('-')) T(1)=R(1) C(NSAMP-1)=N(NSAMP) DO 20 I=2,NSAMP T(I)=T(I-1)+R(I)-M(I) J=NSAMP-I 20 C(J)=C(J+1)+N(J+1)-R(J+1) Z=0.0 K=0 DO 30 I=2,NSAMP-1 P1=.0D0 P2=.0D0 IF(T(I-1).GE..0D0.AND.C(I-1).GE..0D0)THEN P1=M(I)/T(I-1) P2=N(I)/C(I-1) Z1=0.0 DENOM=N(I)*(C(I-1)-T(I-1))*(C(I-1)-N(I))/ / (T(I-1)*C(I-1)**2*(C(I-1)-.1D1)) IF(C(I-1).NE.T(I-1))Z1=(P1-P2)/DSQRT(DENOM) Z=Z+Z1 WRITE(7,25)I,P1,P2,Z1,ZPROB(Z1) 25 FORMAT(1X,I10,4(5X,F10.4)) K=K+1 ENDIF 30 CONTINUE XK=K Z=Z/DSQRT(XK) IF(K.GT.0)WRITE(7,40)Z,ZPROB(Z) 40 FORMAT(/' Overall test statistic : Z = ',F10.4,' Probability = ', , F8.4) IF(K.LE.0)WRITE(7,50) 50 FORMAT(/' Test not performed - insufficeint data!') RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE TDBL1(NTIME,R,S,XD,XB,DBL1,DBL1DF,M,T) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER I,M,N1,NTIME,J,DBL1DF,IDIM PARAMETER (IDIM=50) LOGICAL POOLD DIMENSION R(IDIM),S(IDIM),XD(IDIM),XB(IDIM),A(6),T(6,IDIM) DATA A/6*.0D0/ C C SUBROUTINE COMPUTES Component L1(I) FOR THE TEST BETWEEN C MODELS D AND B (WITH POOLING). C N1=NTIME-1 POOLD=.FALSE. DO 40 I=1,N1 IF(XD(I).LT..0D0.OR.XD(I).GT..1D1)GO TO 40 IF(XB(I).LT..0D0.OR.XB(I).GT..1D1)GO TO 40 A(1)=A(1)+R(I) A(2)=A(2)+S(I) A(3)=A(3)+(.1D1-XD(I))*S(I) A(4)=A(4)+S(I)*XD(I) A(5)=A(5)+(.1D1-XB(I))*S(I) A(6)=A(6)+S(I)*XB(I) IF(A(3).LT..2D1.OR.A(4).LT..2D1.OR. . A(5).LT..2D1.OR.A(6).LT..2D1)GO TO 30 IF(I.NE.NTIME-2)GO TO 10 IF(S(N1)*(.1D1-XD(N1)).LT..2D1.OR.S(N1)*XD(N1).LT..2D1.OR. . S(N1)*(.1D1-XB(N1)).LT..2D1.OR.S(N1)*XB(N1).LT..2D1)GO TO 30 10 DBL1DF=DBL1DF+1 IF(POOLD)T(M,I)=A(2)*(A(6)-A(4))*(A(6)-A(4))/A(4)/A(3) IF(.NOT.POOLD)T(M,I)=S(I)*(XB(I)-XD(I))*(XB(I)-XD(I))/XD(I)/ / (.1D1-XD(I)) POOLD=.FALSE. T(M,IDIM)=T(M,IDIM)+T(M,I) DO 20 J=1,6 20 A(J)=.0D0 GO TO 40 30 POOLD=.TRUE. 40 CONTINUE DBL1DF=NTIME-1-DBL1DF RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE TDBL2(NTIME,M,Z,PD,PB,XD,XB,DBL2,DBL2DF,MD,T) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER I,MD,N1,NTIME,J,DBL2DF,IDIM PARAMETER (IDIM=50) LOGICAL POOLD DIMENSION M(IDIM),Z(IDIM),PB(IDIM),XD(IDIM),XB(IDIM),A(6), , T(6,IDIM) DATA A/6*.0D0/ RHOD(I)=PD /(.1D1-(.1D1-PD )*XD(I)) RHOB(I)=PB(I)/(.1D1-(.1D1-PB(I))*XB(I)) C C SUBROUTINE COMPUTES Component L2(I) FOR THE TEST BETWEEN C MODELS D AND B (WITH POOLING). C POOLD=.FALSE. N1=NTIME-1 DO 40 I=2,NTIME-1 IF(RHOD(I).LT..0D0.OR.RHOD(I).GT..1D1)GO TO 40 IF(RHOB(I).LT..0D0.OR.RHOB(I).GT..1D1)GO TO 40 A(1)=A(1)+M(I) A(2)=A(2)+Z(I) A(3)=A(3)+(M(I)+Z(I))*RHOD(I) A(4)=A(4)+(M(I)+Z(I))*(.1D1-RHOD(I)) A(5)=A(5)+(M(I)+Z(I))*RHOB(I) A(6)=A(6)+(M(I)+Z(I))*(.1D1-RHOB(I)) IF(A(3).LT..2D1.OR.A(4).LT..2D1.OR. . A(5).LT..2D1.OR.A(6).LT..2D1)GO TO 30 IF(I.NE.NTIME-2)GO TO 10 IF(RHOD(N1)*(M(N1)+Z(N1)).LT..2D1.OR. . (.1D1-RHOD(N1))*(M(N1)+Z(N1)).LT..2D1.OR. . RHOB(N1)*(M(N1)+Z(N1)).LT..2D1.OR. . (.1D1-RHOB(N1))*(M(N1)+Z(N1)).LT..2D1)GO TO 30 10 DBL2DF=DBL2DF+1 IF(POOLD)T(MD,I)=(A(1)+A(2))*(A(5)-A(3))**2/A(3)/A(4) IF(.NOT.POOLD)T(MD,I)=(M(I)+Z(I))*(RHOB(I)-RHOD(I))**2/ / RHOD(I)/(.1D1-RHOD(I)) POOLD=.FALSE. T(MD,IDIM)=T(MD,IDIM)+T(MD,I) DO 20 J=1,6 20 A(J)=.0D0 GO TO 40 30 POOLD=.TRUE. 40 CONTINUE DBL2DF=NTIME-2-DBL2DF RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE TDA1(NTIME,R,S,XD,DAL1,DAT1,DAL1DF,MODB,T) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER NTIME,N1,I,J,MODB,DAL1DF,IDIM PARAMETER (IDIM=50) DIMENSION R(IDIM),S(IDIM),XD(IDIM),A(4),T(6,IDIM) C C SUBROUTINE COMPUTES Component L1(I) FOR THE TEST BETWEEN C MODELS D AND A (WITH POOLING). C (ALSO COMPUTES L1(I) FOR TEST OF MODEL B VS A) C N1=NTIME-1 DO 10 I=1,4 10 A(I)=.0D0 DO 40 I=1,NTIME-1 IF(XD(I).LT..0D0.OR.XD(I).GT..1D1)GO TO 40 A(1)=A(1)+R(I) A(2)=A(2)+S(I)-R(I) A(3)=A(3)+(.1D1-XD(I))*S(I) A(4)=A(4)+XD(I)*S(I) IF(A(1).LT..2D1.OR.A(2).LT..2D1.OR. . A(3).LT..2D1.OR.A(4).LT..2D1)GO TO 40 IF(I.NE.NTIME-2)GO TO 20 IF(R(N1).LT..2D1.OR.(S(N1)-R(N1)).LT..2D1.OR. . S(N1)*(.1D1-XD(N1)).LT..2D1.OR.S(N1)*XD(N1).LT..2D1)GO TO 40 20 T(MODB,I)=(A(2)+A(1))*(A(1)-A(3))**2 /A(4)/A(3) DAL1DF=DAL1DF+1 T(MODB,IDIM)=T(MODB,IDIM)+T(MODB,I) DO 30 J=1,4 30 A(J)=.0D0 40 CONTINUE DAL1DF=NTIME-1-DAL1DF RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE TDA2(NTIME,M,Z,PD,XD,DAL2,DAT2,DAL2DF,MODB,T) IMPLICIT DOUBLE PRECISION (A-Z) INTEGER NTIME,N1,I,J,MODB,DAL2DF,IDIM PARAMETER (IDIM=50) DIMENSION M(IDIM),Z(IDIM),PD(1),XD(IDIM),A(4),T(6,IDIM) C C SUBROUTINE COMPUTES Component L2(I) FOR THE TEST BETWEEN C MODELS D AND A (WITH POOLING). C (ALSO COMPUTES L2(I) FOR TEST OF MODEL B VS A) C N1=NTIME-1 DO 10 I=1,4 10 A(I)=.0D0 DO 40 I=2,N1 RHOD=PD(1)/(.1D1-(.1D1-PD(1))*XD(I)) IF(MODB.EQ.4)RHOD=PD(I)/(.1D1-(.1D1-PD(I))*XD(I)) IF(RHOD.LT..0D0.OR.RHOD.GT..1D1)GO TO 40 A(1)=A(1)+M(I) A(2)=A(2)+Z(I) A(3)=A(3)+(M(I)+Z(I))*RHOD A(4)=A(4)+(M(I)+Z(I))*(.1D1-RHOD) IF(A(1).LT..2D1.OR.A(2).LT..2D1.OR. . A(3).LT..2D1.OR.A(4).LT..2D1)GO TO 40 IF(I.NE.NTIME-2)GO TO 20 IF(M(N1).LT..2D1.OR.Z(N1).LT..2D1.OR.RHOD*(M(N1)+Z(N1)).LT. . .2D1.OR.(.1D1-RHOD)*(M(N1)+Z(N1)).LT..2D1)GO TO 40 20 T(MODB,I)=(A(1)+A(2))*(A(1)-A(3))**2 /A(3)/A(4) DAL2DF=DAL2DF+1 T(MODB,IDIM)=T(MODB,IDIM)+T(MODB,I) DO 30 J=1,4 30 A(J)=.0D0 40 CONTINUE DAL2DF=NTIME-2-DAL2DF RETURN END C$CONTROL SEGMENT 'VCB' SUBROUTINE VARCVB(S,M,N,Z,FI,P,T,CHI,V,VCFIP,PRT) C C PRINT VAR-COV MATRICES OF THE ESTIMATES OF: C COV(M(I),M(J)) I=2...S;J=I...S C COV(P(I),M(J)) I=2...S;J=2...S C COV(N(I),N(J)) I=2...S-1;J=I+1...S C COV(B(I),B(J)) I=2...S-2;J=I+1...S-1 C C METHOD: SEE BROWNIE PAPER FOR FORMULAE C IMPLICIT DOUBLE PRECISION (A-Z) INTEGER I,J,K,L,S,PRT,IDIM PARAMETER (IDIM=50) DOUBLE PRECISION T(IDIM),M(IDIM),N(IDIM),Z(IDIM),P(IDIM), , CHI(IDIM),V(IDIM,IDIM),VCFIP(IDIM),VCM(IDIM,IDIM), , VCPM(IDIM,IDIM),VCP(IDIM,IDIM),VCN(IDIM,IDIM),VCB(IDIM,IDIM), , ALFA(IDIM,IDIM),X(IDIM,IDIM),W(IDIM,IDIM),XM(IDIM),XN(IDIM), , XB(IDIM) CHARACTER*8 CHR(3) EQUIVALENCE (XM(1),XN(1),XB(1)) EQUIVALENCE (VCM(1,1),X(1,1),VCN(1,1),VCB(1,1)) DATA CHR/' M ',' N ',' B '/ C C STATEMENT FUNCTIONS C Q(I)=.1D1-P(I) PHI(I)=FI**T(I) XU(I)=(N(I)-M(I))/P(I) C C C INITIALIZE ARRAYS C DO 10 I=1,S XM(I)=(M(I)+Z(I))/(.1D1-CHI(I)+CHI(I)*P(I)) DO 10 J=1,S ALFA(J,I)=.0D0 X(J,I)=.0D0 VCM(J,I)=.0D0 10 VCPM(J,I)=.0D0 C C COMPUTE ARRAY ALFA C ALFA(1,1)=.1D1 DO 30 J=2,S-1 ALFA(J,J)=.1D1 DO 20 K=1,J-1 ALFA(K,J)=.1D1 DO 20 L=K+1,J 20 ALFA(K,J)=ALFA(K,J)*Q(L)*PHI(L) 30 CONTINUE C C COMPUTE ARRAY X C DO 50 J=2,S-1 DO 40 K=2,J 40 X(1,J)=X(1,J)+T(K-1)*PHI(K)/CHI(K)*XM(K)*ALFA(K,J) 50 X(1,J)=X(1,J)*Q(J+1)*CHI(J+1)/FI C DO 60 I=2,S 60 X(I,I-1)=-Q(I)*CHI(I)*XM(I)/(CHI(I-1)*(.1D1-CHI(I-1))) C DO 70 I=2,S-1 DO 70 J=I,S-1 70 X(I,J)=Q(J+1)*CHI(J+1)*(.1D1-PHI(I-1))*XM(I)*PHI(I)*ALFA(I,J)/ / (PHI(I-1)*CHI(I-1)*CHI(I)) C C COMPUTE W=V*X (MATRIX MULT.) C CALL MMULT(V,X,W,S,S,S-1,IDIM,IDIM,IDIM) C C COMPUTE VAR-COV MATRIX OF P(I),M(J) C DO 80 I=2,S DO 80 J=2,S TMP=-T(I-1)*(.1D1-Q(I)*CHI(I))/FI T1=TMP*W(1,J-1)-W(I,J-1)/PHI(I-1) T2=TMP*V(1,1)-V(1,I)/PHI(I-1) T3=TMP*V(1,J)-V(I,J)/PHI(I-1) IF(I.EQ.S)GO TO 80 T1=T1+Q(I)*W(I+1,J-1) T2=T2+Q(I)*V(1,I+1) T3=T3+Q(I)*V(J,I+1) 80 VCPM(I,J)=.1D1/CHI(I)*(T1+T(J-1)*XM(J)/FI*T2+ + XM(J)/(.1D1-CHI(J-1))*T3) IF(PRT.EQ.1)WRITE(7,90) 90 FORMAT(/' Variance-covariance matrix: COV[ p(i), M(j) ]', , ' ; I=2,...S ; J=2,...S') DO 100 I=2,S 100 IF(PRT.EQ.1)WRITE(7,110)(VCPM(I,J),J=2,S) 110 FORMAT(/(10(1X,F12.2))) C C COMPUTE VAR-COV MATRIX OF P C IF(PRT.EQ.1)WRITE(7,120) 120 FORMAT(/' Variance-covariance matrix: COV[ p(i), p(j) ]', , ' ; i=2,...s ; j=2,...s') DO 140 I=2,S DO 130 J=I,S T1=T(I-1)*(.1D1-Q(I)*CHI(I))/FI T2=T(J-1)*(.1D1-Q(J)*CHI(J))/FI VCP(I,J)=T1*(T2*V(1,1)+V(1,J)/PHI(J-1))+ + (T2*V(1,I)+V(I,J)/PHI(J-1))/PHI(I-1) IF(I.LT.S)VCP(I,J)=VCP(I,J)- - Q(I)*(T2*V(1,I+1)+V(I+1,J)/PHI(J-1)) IF(J.LT.S)VCP(I,J)=VCP(I,J)- - Q(J)*(T1*V(1,J+1)+V(I,J+1)/PHI(I-1)) IF(I.LT.S.AND.J.LT.S)VCP(I,J)=VCP(I,J)+ + Q(I)*Q(J)*V(I+1,J+1) VCP(I,J)=VCP(I,J)/(CHI(I)*CHI(J)) 130 VCP(J,I)=VCP(I,J) 140 IF(PRT.EQ.1)WRITE(7,240)(VCP(I,J),J=2,S) C C COMPUTE VAR-COV MATRIX OF M C IF(PRT.EQ.1)WRITE(7,150) 150 FORMAT(/' Variance-covariance matrix: COV[ M(i), M(j) ]', , ' ; i=2,...s; j=2,...s') DO 170 I=2,S DO 160 J=I,S VCM(I,J)=Q(J)*CHI(J)*XM(I)*ALFA(I-1,J-1)/ / (.1D1-Q(I)*CHI(I))+ + XM(I)*T(I-1)/FI*W(1,J-1)+ + XM(I)/(.1D1-CHI(I-1))*W(I,J-1)+ + XM(J)*T(J-1)/FI* * (W(1,I-1)+XM(I)*T(I-1)/FI*V(1,1)+ + XM(I)/(.1D1-CHI(I-1))*V(1,I))+ + XM(J)/(.1D1-CHI(J-1))* * (W(J,I-1)+XM(I)*T(I-1)/FI*V(1,J)+ + XM(I)/(.1D1-CHI(I-1))*V(J,I)) 160 VCM(J,I)=VCM(I,J) 170 IF(PRT.EQ.1)WRITE(7,110)(VCM(I,J),J=2,S) C CALL OUTVAR(XM,2,S,VCM,CHR(1)) C C COMPUTE VAR-COV MATRIX OF N C IF(PRT.EQ.1)WRITE(7,180) 180 FORMAT(/' Variance-covariance matrix: COV[ N(i), N(j) ]', , ' ; i=2,...s; j=2,...s') DO 200 I=2,S XN(I)=XU(I)+XM(I) DO 190 J=I,S VCN(I,J)=VCM(I,J)-XU(I)/P(I)*VCPM(I,J)- - XU(J)/P(J)*VCPM(J,I)+ + XU(I)*XU(J)/(P(I)*P(J))*VCP(I,J) 190 VCN(J,I)=VCN(I,J) 200 IF(PRT.EQ.1)WRITE(7,110)(VCN(I,J),J=2,S) C CALL OUTVAR(XN,2,S,VCN,CHR(2)) C C COMPUTE VAR-COV MATRIX: COV[B(I),B(J)] C IF(PRT.EQ.1)WRITE(7,210) 210 FORMAT(/' Variance-covariance matrix: COV[ B(i), B(j) ]', , ' ; i=2,...s-1; j=2,...s-1') DO 230 I=2,S-1 DO 220 J=I,S-1 VCB(I,J)=Q(I)*PHI(I)*T(I)*XU(I)/FI* * (Q(J)*PHI(J)*T(J)*XU(J)/FI*V(1,1)-PHI(J)*XU(J)/P(J)*VCFIP(J)+ + XU(J+1)/P(J+1)*VCFIP(J+1))+PHI(I)*XU(I)/P(I)* * (-Q(J)*PHI(J)*T(J)*XU(J)/FI*VCFIP(I)+PHI(J)*XU(J)/P(J)*VCP(I,J)- - XU(J+1)/P(J+1)*VCP(I,J+1))+XU(I+1)/P(I+1)* * (Q(J)*PHI(J)*T(J)*XU(J)/FI*VCFIP(I+1)- - PHI(J)*XU(J)/P(J)*VCP(I+1,J)+XU(J+1)/P(J+1)*VCP(I+1,J+1)) IF(J.EQ.I+1)VCB(I,J)=VCB(I,J)-Q(J)*PHI(J)*XU(J)/P(J) IF(I.EQ.J)VCB(I,I)=VCB(I,I)+XU(I+1)*Q(I+1)/P(I+1)+ + Q(I)*PHI(I)*XU(I)*(.1D1+Q(I)*PHI(I)/P(I)) VCB(J,I)=VCB(I,J) 220 CONTINUE IF(PRT.EQ.1)WRITE(7,110)(VCB(I,J),J=2,S-1) 230 XB(I)=XU(I+1)-PHI(I)*Q(I)*XU(I) C CALL OUTVAR(XB,2,S-1,VCB,CHR(3)) 240 FORMAT(/(10(1X,F12.8))) RETURN END C$CONTROL SEGMENT 'MAIN' SUBROUTINE OUTVAR(VAR,I1,I2,VCMAT,CHR) PARAMETER (IDIM=50) DOUBLE PRECISION DSQRT,C1,SE,C2,TITLE,VAR(IDIM),VCMAT(IDIM,IDIM) DOUBLE PRECISION MEAN,SEMEAN CHARACTER*8 CHR INTEGER S,YR COMMON /BLK0/S,YR(IDIM),TITLE(16) MEAN=.0D0 SEMEAN=.0D0 C C OUTPUT VARIABLE WITH 95% CONFIDENCE INTERVAL C WRITE(7,10)CHR C 10 FORMAT(/38X,'Standard'/' Period',7X,A8,5X,'Variance',4X,'error', ,8X,' 95% confidence interval'/' ',75('-')) DO 30 I=I1,I2 MEAN=MEAN+VAR(I) DO 20 J=I1,I2 20 SEMEAN=SEMEAN+VCMAT(I,J) SE=.1D30 IF(VCMAT(I,I).GE..0D0)SE=DSQRT(VCMAT(I,I)) C1=VAR(I)-SE*.196D1 C2=VAR(I)+SE*.196D1 WRITE(7,40)YR(I),VAR(I),VCMAT(I,I),SE,C1,C2 30 IF(MOD(YR(I),5).EQ.0)WRITE(7,40) 40 FORMAT(I5,4X,F11.2,3X,F12.4,2X,F8.2,6X,F11.2,' - ',F11.2) MEAN=MEAN/FLOAT(I2-I1+1) SEMEAN=SEMEAN/FLOAT(I2-I1+1) SE=.1D30 IF(SEMEAN.GT..0D0)SE=DSQRT(SEMEAN) C1=MEAN-SE*.196D1 C2=MEAN+SE*.196D1 WRITE(7,50)MEAN,SEMEAN,SE,C1,C2 50 FORMAT(1X,75('-')/' MEAN',4X,F11.2,3X,F12.4,2X,F8.2,6X,F11.2, , ' - ',F11.2/1X,75('-')) RETURN END C$CONTROL SEGMENT 'VCD' SUBROUTINE VARCVD(FI,P,M,N,NS,Z,T,X,V,PRT) C C C COMPUTE & PRINT VAR-COV MATRICES FOR MODEL D C IMPLICIT DOUBLE PRECISION (A-Z) INTEGER I,J,K,L,S,YR,PRT,IDIM PARAMETER (IDIM=50) DIMENSION T(IDIM),M(IDIM),N(IDIM),Z(IDIM),X(IDIM),XX(2,IDIM), , W(IDIM+1),WW(2,IDIM),V(2,2),ALFA(IDIM,IDIM),SS(IDIM),PP(IDIM), , XM(IDIM),XN(IDIM),U(IDIM),NS(IDIM),COVM(IDIM,IDIM),COVPM(IDIM), , COVN(IDIM,IDIM),B(IDIM),COVB(IDIM,IDIM) CHARACTER*8 CHR(3) EQUIVALENCE (COVM(1,1),COVN(1,1),COVB(1,1)) COMMON /BLK0/S,YR(IDIM),TITLE(16) DATA CHR/' M ',' N ',' B '/ C C STATEMENT FUNCTION PHI(I)=FI**T(I) C C COMPUTE W,ALFA C Q=.1D1-P W(1)=.0D0 DO 30 J=1,S U(J)=(N(J)-M(J))/P XM(J)=(M(J)+Z(J))/(.1D1-Q*X(J)) IF(J.NE.S)W(J+1)=PHI(J)**2*(Q*Q*W(J)+NS(J)/X(J)) DO 30 K=1,J SUM=.0D0 IF(K.EQ.J)GO TO 20 DO 10 L=K+1,J 10 SUM=SUM+T(L) 20 ALFA(K,J)=Q**(J-K)*FI**SUM 30 CONTINUE C C COMPUTE 2ND X: XX C XX(2,1)=.0D0 DO 100 J=2,S SUM1=.0D0 SUM2=.0D0 DO 40 K=2,J 40 SUM1=SUM1+T(K-1)*ALFA(K-1,J-1)*(XM(K)+(.1D1-Q*X(K))*W(K)) IF(J.EQ.S)GO TO 60 DO 50 K=J+1,S 50 SUM2=SUM2+T(K-1)*ALFA(J-1,K-1)*(.1D1-Q*X(K)) 60 XX(1,J-1)=Q*X(J)/FI*(SUM1+(XM(J)/(.1D1-Q*X(J))+W(J))*SUM2) SUM1=.0D0 SUM2=.0D0 IF(J.LE.2)GO TO 80 DO 70 K=2,J-1 70 SUM1=SUM1+ALFA(K,J-1)*PHI(K)*(Q*X(K)*W(K)-XM(K)) 80 DO 90 K=J,S 90 SUM2=SUM2+X(K)*ALFA(J-1,K-1) 100 XX(2,J-1)=Q*X(J)*SUM1+Q*X(J)*(XM(J)/(.1D1-Q*X(J))+W(J))*SUM2 C C COMPUTE 2ND W:WW = V * X C CALL MMULT(V,XX,WW,2,2,S-1,2,2,IDIM) C C COMPUTE 2ND S,P:SS,PP & VAR-COV MATRIX FOR M C IF(PRT.EQ.1)WRITE(7,120) 120 FORMAT(/' Variance-covariance matrix: COV[ M(i), M(j) ]', , ' ; i=2,...s; j=2,...s') SS(S)=.0D0 PP(S)=.1D1 DO 140 I=2,S-1 SUM1=.0D0 SUM2=.0D0 DO 130 K=I,S-1 SUM1=SUM1+T(K)*ALFA(I-1,K-1)*(.1D1-X(K)) 130 SUM2=SUM2+ALFA(I,K)*X(K+1) SS(I)=Q/FI*SUM1 140 PP(I)=X(I)+Q*PHI(I)*SUM2 DO 170 I=2,S DO 150 J=I,S QI=.1D1-Q*X(I) QJ=.1D1-Q*X(J) COVM(I,J)=Q*X(J)*XM(I)/QI*ALFA(I-1,J-1)- - XM(J)/QJ*(SS(J)*WW(1,I-1)+PP(J)*WW(2,I-1))- - XM(I)*SS(I)/QI*(WW(1,J-1)-SS(J)*XM(J)/QJ*V(1,1)- - PP(J)*XM(J)/QJ*V(1,2))- - XM(I)*PP(I)/QI*(WW(2,J-1)-SS(J)*XM(J)/QJ*V(1,2)- - PP(J)*XM(J)/QJ*V(2,2)) 150 COVM(J,I)=COVM(I,J) IF(PRT.EQ.1)WRITE(7,160)(COVM(I,J),J=2,S) 160 FORMAT(10F13.2) 170 CONTINUE C CALL OUTVAR(XM,2,S,COVM,CHR(1)) C C COMPUTE COV(P,M(I)) C DO 180 I=2,S QI=.1D1-Q*X(I) 180 COVPM(I)=WW(2,I-1)-XM(I)*SS(I)/QI*V(1,2)-XM(I)*PP(I)/QI*V(2,2) IF(PRT.EQ.1)WRITE(7,190) 190 FORMAT(/' Variance-covariance matrix: COV[ N(i), N(j) ]', , ' ; i=2,...s; j=2,...s') DO 210 I=2,S XN(I)=U(I)+XM(I) DO 200 J=I,S COVN(I,J)=COVM(I,J)-U(J)/P*COVPM(I)-U(I)/P*COVPM(J)+ + U(I)*U(J)/P**2*V(2,2) 200 COVN(J,I)=COVN(I,J) COVN(I,I)=COVN(I,I)+U(I)/P*Q 210 IF(PRT.EQ.1)WRITE(7,160)(COVN(I,J),J=2,S) C CALL OUTVAR(XN,2,S,COVN,CHR(2)) C C COMPUTE VAR-COV (B) C IF(PRT.EQ.1)WRITE(7,220) 220 FORMAT(/' Variance-covariance matrix: COV[ B(i), B(j) ]', , ' ; i=2,...s-1; j=2,...s-1') DO 240 I=2,S-1 B(I)=U(I+1)-PHI(I)*Q*U(I) DO 230 J=I,S-1 COVB(I,J)=V(2,2)*((U(I+1)-PHI(I)*U(I))/P* * (U(J+1)-PHI(J)*U(J))/P)+ + V(1,1)*(T(I)*Q*PHI(I)*U(I)/FI* * T(J)*Q*PHI(J)*U(J)/FI)+ + V(1,2)*(T(I)*Q*PHI(I)*U(I)/FI* * (U(J+1)-PHI(J)*U(J))/P+ + T(J)*Q*PHI(J)*U(J)/FI* * (U(I+1)-PHI(I)*U(I))/P) IF(J.EQ.I+1)COVB(I,J)=COVB(I,J)-PHI(I+1)*Q*U(I+1)/P 230 COVB(J,I)=COVB(I,J) COVB(I,I)=COVB(I,I)+Q*U(I+1)/P+Q*PHI(I)*U(I)*(.1D1+Q*PHI(I)/P) 240 IF(PRT.EQ.1)WRITE(7,160)(COVB(I,J),J=2,S-1) C CALL OUTVAR(B,2,S-1,COVB,CHR(3)) C RETURN END