SUBROUTINE SIMPLE(K,T,MATRIX,N,YRONE,HEADER)
C *** SIMPLE COMPUTES ESTIMATES OF THE CONSTANT SURVIVAL AND RECOVERY
C *** RATES USING THE ITERATIVE METHODS EMPLOYED IN MODEL 3.  SIMPLE
C *** IS USEFUL (AS AN OPTION--SEE MAIN) IF THE DATA ARE VERY POOR, E.G.
C *** NO BANDINGS IN SOME YEARS, NO RECOVERIES FROM SMALL BANDED
C *** SAMPLES, ETC.  NO TESTING IS DONE AND THE VARIANCES ARE USUALLY
C *** SUBSTANTIALLY UNDERESTIMATED.
      IMPLICIT DOUBLE PRECISION(A-G,O-Z)
      CHARACTER*80 HEADER
      DOUBLE PRECISION INFO,NEXT,HIP,HIL,LOWP,LOWL,INVERS,LAMBDA,N,
     1 MATRIX,LOWF,HIF
      DOUBLE PRECISION MLS,LOWMLS,HIMLS
      INTEGER YEAR,YRONE,BANDS,T,IMAT
      DIMENSION MATRIX(21,21),N(20)
      COMMON /SIMPCM/ R(20),A(20),PD(2),AHAT(20),CHISQS(21
     1,21),INVERS(2,2),INFO(2,2),NEXT(2),THCHI(100),EE(21,21),IMAT(20),
     2 EREC(20),ISAVE(20)
      DO 5 I=1,T
      R(I)=0.0
      DO 5 II=1,K
      R(I)=R(I)+MATRIX(I,II)
    5 CONTINUE
      ITER=0
      BIGN=0.0
      BIGT=0.0
      Q=0.0
      DO 10 I=1,T
      BIGT=BIGT+R(I)
      BIGN=BIGN+N(I)
      A(I)=N(I)
   10 CONTINUE
      DO 15 I=1,T
      DO 15 J=1,K
      Q=Q+(J-I)*MATRIX(I,J)
   15 CONTINUE
      LAMBDA=BIGT/BIGN
      PHI=Q/(BIGT+Q)
      NEXT(1)=LAMBDA
      NEXT(2)=PHI
      WRITE (6,20)
   20 FORMAT ('0','ANALYSIS ASSUMING CONSTANT SURVIVAL AND RECOVERY ',
     1'PROBABILITIES')
      WRITE (6,25)
   25 FORMAT (1X,'(A GENERALIZATION AND EXTENSION OF THE MODELS DEVELOP'
     1,'ED BY CHAPMAN AND ROBSON (1960. BIOMETRICS) AND ')
      WRITE (6,30)
   30 FORMAT (1X,' HALDANE (1955. PROC. XI INT. ORN. CONGR.) -- SEE BOT'
     1,'TOM OF PAGE 245 OF BOOK BY SEBER))'/)
      WRITE (6,35) HEADER
   35 FORMAT (/' ',20X,A//)
      WRITE (6,40)
   40 FORMAT (' ','YEAR NUMBER',35X,'RECOVERY MATRIX')
      WRITE (6,45)
   45 FORMAT (' ',5X,'BANDED')
      WRITE (6,50)
   50 FORMAT ('+','____ ______',35X,'_______________'/)
      YEAR=YRONE-1
      DO 65 I=1,T
      BANDS=A(I)
      YEAR=YEAR+1
      DO 55 JJ=1,K
      IMAT(JJ)=MATRIX(I,JJ)
   55 CONTINUE
      WRITE (6,60) YEAR,BANDS,(IMAT(JJ),JJ=1,K)
   60 FORMAT (' ',I4,1X,I5,2X,20(2X,I4))
   65 CONTINUE
      WRITE (6,70)
   70 FORMAT (/' ','INTERMEDIATE STATISTICS')
      WRITE (6,75) BIGN
   75 FORMAT (' ','N  = ',F6.0)
      WRITE (6,80) BIGT
   80 FORMAT (' ','T = ',F6.0)
      WRITE (6,85) Q
   85 FORMAT (' ','Q  = ',F6.0/)
   90 CONTINUE
      IF (NEXT(2).GT.1.0) GO TO 145
      IF (NEXT(2).LT.0.) GO TO 145
      IF (NEXT(1).LT.0.) GO TO 145
      ITER=ITER+1
      IF (ITER.GT.40) GO TO 145
      EREC(1)=(1.0-PHI)*LAMBDA
      DO 16 I=2,K
   16 EREC(I)=EREC(I-1)*PHI
      EXPT=0.
      EXPQ=0.
      DO 17 I=1,T
       DO 17 J=I,K
      EXPT=EXPT+A(I)*EREC(J-I+1)
   17  EXPQ=EXPQ+A(I)*(J-I)*EREC(J-I+1)
      SUM1=0.0
      SUM2=0.0
      SUM3=0.0
      SUM4=0.0
      SUM5=0.0
      DO 95 I=1,T
      SUM1=SUM1+(N(I)-R(I))/(1.-LAMBDA*(1.-PHI**(K-I+1)))
      SUM2=SUM2+(N(I)-R(I))*(K-I+1)*(PHI**(K-I))/(1.-LAMBDA*(1.-PHI**(K-
     1I+1)))
      SUM3=SUM3+(N(I)*(K-I+1)*PHI**(K-I))/(1.-LAMBDA*(1.-PHI**(K-I+1)))
      SUM4=SUM4+(N(I)/(1.-LAMBDA*(1.-PHI**(K-I+1))))
      SUM5=SUM5-LAMBDA*N(I)*(K-I+1)*(K-I)*(PHI**(K-I-1))
      AA=N(I)*LAMBDA*LAMBDA*(K-I+1)*(K-I+1)
      AA=AA*(PHI**(2.*(K-I)))
      AA=AA/(1.-LAMBDA*(1.-PHI**(K-I+1)))
      SUM5=SUM5+AA
   95 CONTINUE
      SUM3=-1*SUM3
      SUM4=SUM4-BIGN
      PD(1)=(BIGN/LAMBDA)-(1./LAMBDA)*SUM1
      PD(2)=((-1*BIGT)/(1.-PHI))+(Q/PHI)+(LAMBDA*SUM2)
      INFO(1,2)=SUM3
      INFO(2,1)=SUM3
      INFO(1,1)=((1./LAMBDA)**2)*SUM4
      INFO(2,2)=(EXPT/((1.-PHI)**2))+(EXPQ/PHI**2)+SUM5
C     COMPUTE THE INVERSE OF THE INFORMATION MATRIX
      DET=INFO(1,1)*INFO(2,2)-(INFO(2,1)*INFO(1,2))
      INVERS(1,1)=INFO(2,2)/DET
      INVERS(2,2)=INFO(1,1)/DET
      INVERS(2,1)=(-1.*INFO(2,1))/DET
      INVERS(1,2)=(-1.*INFO(1,2))/DET
C     THE INTERATION IS:
      NEXT(1)=LAMBDA+((INVERS(1,1)*PD(1))+(INVERS(1,2)*PD(2)))
      NEXT(2)=PHI+((INVERS(2,1)*PD(1))+(INVERS(2,2)*PD(2)))
      IF (DABS(NEXT(1)-LAMBDA).LT.0.00001.AND.DABS(NEXT(2)-PHI).LT.
     +0.00001) GO TO 100
      LAMBDA=NEXT(1)
      PHI=NEXT(2)
      GO TO 90
  100 NEXT(1)=NEXT(1)*100.
      NEXT(2)=NEXT(2)*100.
      IF (INVERS(2,2).GT.0.0) STDEVP=DSQRT(INVERS(2,2))*100.
      IF (INVERS(1,1).GT.0.0) STDEVL=DSQRT(INVERS(1,1))*100.
      CVP=STDEVP/NEXT(2)*100.
      CVL=STDEVL/NEXT(1)*100.
      LOWP=(NEXT(2)-(1.96*STDEVP))
      HIP=(NEXT(2)+(1.96*STDEVP))
      LOWL=(NEXT(1)-(1.96*STDEVL))
      HIL=(NEXT(1)+(1.96*STDEVL))
      F=(100.-NEXT(2))*(NEXT(1)/100.)
      F=F/100.
      COVAR=INVERS(2,1)
      AMORT=1.-(NEXT(2)/100.)
      VARF=(F/AMORT)**2*INVERS(2,2)-2.*F*COVAR+AMORT**2*INVERS(1,1)
      COVFS=(-1.*(F/AMORT))*INVERS(2,2)+COVAR*AMORT
      F=F*100.
      IF (VARF.GT.0.0) STDEVF=DSQRT(VARF)*100.
      LOWF=F-(1.96*STDEVF)
      HIF=F+(1.96*STDEVF)
      IF (VARF*INVERS(2,2).GT.0.0) CORRFS=COVFS/DSQRT(VARF*INVERS(2,2))
      CVF=(STDEVF/F)*100.
      WRITE (6,105)
  105 FORMAT (/' ',1X,'PARAMETER',31X,'ESTIMATE (%)',07X,'STD. DEV. (%)'
     1,06X,'COEF. VARIAT. (%)',05X,'90% CONFIDENCE INTERVAL')
      WRITE (6,110)
  110 FORMAT (' ',124('-'))
      WRITE (6,115) NEXT(2),STDEVP,CVP,LOWP,HIP
  115 FORMAT (/' ','SURVIVAL PROBABILITY  (S) ',18X,F5.2,13X,F5.2,15X,F5
     1.2,15X,F5.2,' -- ',F5.2/)
      WRITE (6,120) F,STDEVF,CVF,LOWF,HIF
  120 FORMAT (1X,'RECOVERY PROBABILITY (F)             ',7X,F5.2,13X,
     1 F5.2,15X,F5.2,15X,F5.2,' -- ',F5.2)
      WRITE (6,110)
      WRITE (6,125)
  125 FORMAT (1X,'THE STANDARD DEVIATIONS FOR S AND F ARE UNDERESTIMATE'
     1,'D BY A FACTOR OF PERHAPS 1.5 TO 5.0 ')
      CORR=(INVERS(2,1))/((STDEVP*STDEVL)/10000.)
      WRITE (6,130) CORRFS
  130 FORMAT (' ','CORRELATION(S,F) = ',F10.8)
      MLS=-1.0/DLOG(NEXT(2)/100.0)
      SDMLS=DSQRT((1.0/MLS**2)*(1.0/(DLOG(MLS)**4))*INVERS(2,2))
      LOWMLS=MLS-1.96*SDMLS
      HIMLS=MLS+1.96*SDMLS
      WRITE (6,135) MLS
  135 FORMAT (' ','MEAN LIFE SPAN AS AN ADULT = ',F7.2)
      WRITE (6,140) ITER
  140 FORMAT (' ',I5,'  ITERATIONS')
      RETURN
  145 WRITE (6,150)
  150 FORMAT (' ','ESTIMATION NOT POSSIBLE--CHECK THE DATA AND RESUBMIT'
     1)
      RETURN
      END
      SUBROUTINE CHI1(K,Z)
C *** CHI1 AND CHI2 COMPUTE THE PROBABILITY LEVEL OF A CHI SQUARE VALUE
C *** GIVEN ITS DEGREES OF FREEDOM
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON/BLOCK/P,LITE
      DATA PI/1.7724539D0/
C     K=DEGREES OF FREEDOM
C     Z=CHI-SQUARE VALUE
C     P=PROBABILITY INTEGRAL TRANSFORM
C     IT=0 IF THE SUBROUTINE FAILED
      IF (Z.GT.600.0) P=1.0
      IF (Z.GT.600.0) GO TO 35
      P=0.
      IT=0
      IF (Z.LE.0.) RETURN
      IF (K.LE.0) RETURN
      IF (K.GT.100) GO TO 25
      RK=K*.5
      Z2=Z*.5
C
      IRK=RK
      IT=1
      IF (IRK.NE.RK) GO TO 10
C
      Q=DEXP(-Z2)
      P=1.-Q
      IF (IRK.EQ.1) GO TO 35
      MAX=IRK-1
      DO 5 I=1,MAX
      X=I
      Q=Q*Z2/X
    5 P=P-Q
      GO TO 35
C
   10 P=1.
      IF (Z.GT.15.) GO TO 15
C
      CALL CHI2(Z2,P)
   15 IF (K.EQ.1) GO TO 35
C
      A=DSQRT(Z2)*2.*DEXP(-Z2)/PI
      P=P-A
      IF (K.EQ.3) GO TO 35
      I=IRK-2
   20 A=A*Z2/(RK-I-1)
      P=P-A
      I=I-1
      IF (I.LT.0) GO TO 35
      GO TO 20
C
   25 Y=DSQRT(2.*Z)-SQRT(2.*K)
      IT=1
      SIGN=-1.
      IF (Y.GE.0.) SIGN=+1.
      W=Y*Y
      Z2=W*.5
      P=1.
      IF (W.GT.15.) GO TO 30
      CALL CHI2(Z2,P)
   30 P=.5+SIGN*P*.5
      GO TO 35
C
   35 P=1.-P
      IF (IT.NE.0.AND.LITE.NE.1) WRITE (6,40) Z,P
   40 FORMAT (/' ','PROBABILITY OF A CHI-SQUARE VALUE LARGER THAN',F8.2,
     12X,'=',F13.8//)
      RETURN
      END
      SUBROUTINE CHI2(Z2,P)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/1.7724539D0/
      A=DSQRT(Z2)/PI
      I=0
      S=2.
      Q=0.
      FAC=1.
    5 I=I+1
      FAC=-FAC*Z2/I
      S=S+FAC/(I+.5)
      P=A*S
      IF (DABS(P-Q).LT..000000001) RETURN
      Q=P
      GO TO 5
      END
      SUBROUTINE TEST(N,MATRIX,A,THCHI)
C *** TEST IS A SIMPLE TEST OF THE NULL HYPOTHESIS THAT THE FIRST-YEAR
C *** RECOVERY RATES (UNDER MODELS 0 AND 1) ARE CONSTANT.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z,M)
      INTEGER DF
      DIMENSION A(20)  ,TR(20),TC(20),MATRIX(21,21),THCHI(100),B(20,2)
      DO 5 I=1,N
      B(I,1)=MATRIX(I,I)
      B(I,2)=A(I)-MATRIX(I,I)
    5 CONTINUE
      CS=0.0
      DF=N-1
C     ROW TOTALS
      DO 10 I=1,N
      TR(I)=0.0
      DO 10 J=1,2
      TR(I)=TR(I)+B(I,J)
   10 CONTINUE
C     COLUMN TOTALS
      DO 15 J=1,2
      TC(J)=0.0
      DO 15 I=1,N
      TC(J)=TC(J)+B(I,J)
   15 CONTINUE
      GT=0.0
      DO 20 I=1,N
   20 GT=GT+TR(I)
      DO 25 J=1,2
      DO 25 I=1,N
      E=TR(I)*TC(J)/GT
   25 CS=CS+(((B(I,J)-E)**2)/E)
      WRITE (6,30)
   30 FORMAT (///1X,'TEST OF THE NULL HYPOTHESIS THAT THE FIRST-YEAR (D'
     1,'IRECT) RECOVERY RATES ARE CONSTANT EACH YEAR:'/)
      WRITE (6,35) CS
   35 FORMAT (' ','   CHI-SQUARED (SAMPLE) = ',F8.2)
      WRITE (6,40) THCHI(DF)
   40 FORMAT (' ','   THEORETICAL CHI-SQUARE VALUE AT THE 5% LEVEL =',F8
     1.2)
      WRITE (6,45) DF
   45 FORMAT (' ','   DEGREES OF FREEDOM = ',I6)
      CALL CHI1(DF,CS)
      RETURN
      END
      SUBROUTINE FANDFP(S,KK,R,C,T,MATRIX,A,YRONE)
C *** FANDFP IS A CONTINGENCY TABLE TEST OF MODEL 0 VS. MODEL 1
      DOUBLE PRECISION MATRIX,R,C,T,OBS,EXPECT,RTOT,CTOT,A,DRATE,AHAT,
     1SIGN,PART,SUM,NORMAL,P
      DOUBLE PRECISION SUMTHR,SUMTWO
      INTEGER S,SM1 ,YRONE
      COMMON/BLOCK/P,LITE
      DIMENSION MATRIX(21,21),OBS(2,2),EXPECT(2,2),RTOT(2),CTOT(2)
      DIMENSION A(20),DRATE(20),AHAT(20),SIGN(20),PART(20)
      DIMENSION R(20),C(20),T(20)
      SUM=0.0
      LITE=1
      MDIFF=KK-S
      IF (MDIFF.EQ.0) SM1=S-1
      IF (MDIFF.GT.0) SM1=S
      IYEAR=YRONE
      PART(1)=0.0
      CHISUM=0.0
      SUMTWO=0.0
C
      DO 5 I=2,SM1
      SIGN(I)=-1.0

      DRATE(I)=MATRIX(I,I)/A(I)
      AHAT(I)=((C(I)-MATRIX(I,I))*(R(I)-MATRIX(I,I)))/(A(I)*(T(I)-R(I)-C
     1(I)+MATRIX(I,I)))
      IF (DRATE(I).GT.AHAT(I)) SIGN(I)=1.0
    5 CONTINUE
      WRITE (6,10)
   10 FORMAT (///1X,'TEST OF THE NULL HYPOTHESIS THAT THE FIRST-YEAR RE'
     1,'COVERY RATES AND/OR SURVIVAL RATES ARE THE SAME AS THOSE FROM')
      WRITE (6,15)
   15 FORMAT (' ','COHORTS BANDED IN PREVIOUS YEARS'/)
      WRITE (6,20)
   20 FORMAT (1X,'THIS IS A TEST OF MODEL 1 (THE NULL HYPOTHESIS) VS. M'
     1,'ODEL 0 (AN ALTERNATIVE HYPOTHESIS)'/)
C
      WRITE (6,25)
   25 FORMAT (' ','       I      CHI-SQUARE      NORMAL(0,1) ')
      DO 40 I=2,SM1
      CHISUM=0.0
      OBS(1,1)=C(I)-MATRIX(I,I)
      OBS(1,2)=T(I)-C(I)-R(I)+MATRIX(I,I)
      OBS(2,1)=MATRIX(I,I)
      OBS(2,2)=R(I)-MATRIX(I,I)
C
      RTOT(1)=OBS(1,1)+OBS(1,2)
      RTOT(2)=OBS(2,1)+OBS(2,2)
      CTOT(1)=OBS(1,1)+OBS(2,1)
      CTOT(2)=OBS(1,2)+OBS(2,2)
      GRAND=RTOT(1)+RTOT(2)
C
      EXPECT(1,1)=RTOT(1)*CTOT(1)/GRAND
      EXPECT(1,2)=RTOT(1)*CTOT(2)/GRAND
      EXPECT(2,1)=RTOT(2)*CTOT(1)/GRAND
      EXPECT(2,2)=RTOT(2)*CTOT(2)/GRAND
C
      DO 30 II=1,2
      DO 30 JJ=1,2
      CHI=((OBS(II,JJ)-EXPECT(II,JJ))**2)/EXPECT(II,JJ)
      CHISUM=CHISUM+CHI
   30 CONTINUE
      SUMTWO=CHISUM+SUMTWO
      PART(I)=SQRT(CHISUM)*SIGN(I)
      IYEAR=IYEAR+1
      WRITE (6,35) IYEAR,CHISUM,PART(I)
   35 FORMAT (' ',4X,I4,5X,F8.2,9X,F8.2)
      SUM=SUM+PART(I)
   40 CONTINUE
C
      WRITE (6,45) SUMTWO,SUM
   45 FORMAT (' ',4X,'TOTAL',4X,F8.2,9X,F8.2/)
      SM1=SM1-1
      NORMAL=SUM**2/SM1
      CALL CHI1(1,NORMAL)
      SIGN1=1.0
      IF (SUM.LT.0.) SIGN1=-1.0
      PROB=(1.0+SIGN1*(1.-P))/2.0
      PROB=1.-PROB
      WRITE (6,50) SUM,PROB
   50 FORMAT (' ','PROBABILITY OF OBSERVING A VALUE LARGER THAN',F8.2,'
     1   IS ',F13.8)
      LITE=0
      SUMTHR=SUMTWO
      CALL CHI1(SM1,SUMTHR)
      RETURN
      END