C------------------------------------------------------------------------
C
C  Program name:  CONEST
C
C  Description:   For integer-valued censored Normal(mu,sigma) data,
C                 estimate mu and sigma from a contour plot of the
C                 log likelihood and plot confidence bands for the
C                 estimates.
C
C  Last revised:  July 1989  (currently under development)
C
C-----------------------------------------------------------------------
C
      REAL*8  MUEST, SIGEST, LMAX
      INTEGER MAXI, NI(10), OPT
C
      LMAX=-999999999.9D0
C
      WRITE(*,10)
   10 FORMAT(////,20X,'CONEST',////)
C
      CALL INPUT(MAXI,NI)
C
      CALL CONPLT(MAXI,NI,MUEST,SIGEST,LMAX)
C
   20 WRITE(*,90)
   90 FORMAT(12(/),1X,'CHOOSE OPTION:  ',
     +////,5X,' (1) PRODUCE PLOT OF LIKELIHOOD ',
     +'OVER ANOTHER REGION',
     +//,5X,' (2) VIEW CONFIDENCE BANDS AROUND BEST ',
     +'ESTIMATES FOUND THIS SESSION',
     +//,5X,' (3) EXIT FROM THE PROGRAM',
     +8(/),1X,'ENTER NUMBER OF SELECTION:  ')
C
      READ(*,*) OPT
      WRITE(*,32)
   32 FORMAT(20(/))
      IF (OPT.GE.1.AND.OPT.LE.3) GOTO 107
C
      WRITE(*,105)
  105 FORMAT(//,1X,'INVALID SELECTION.  PLEASE TRY AGAIN.')
      GOTO 20
C
  107 GOTO (110,120,130),OPT
C
  110 CALL CONPLT(MAXI,NI,MUEST,SIGEST,LMAX)
      GOTO 20
C
  120 CALL BNDPLT(MAXI,NI,MUEST,SIGEST,LMAX)
      GOTO 20
C
  130 STOP
      END
C
C------------------------------------------------------------------
C
      SUBROUTINE INPUT(MAXI,NI)
C
      INTEGER MAXI,NI(10)
C
      WRITE(*,20)
   20 FORMAT(//,1X,'Data values must be coded to integer values ',
     +             'beginning with 1.',
     +      ///,1X,'Enter the maximum integer value in the coded ',
     +             'data and return',/,1X,'---->  ')
      READ(*,*) MAXI
C
      DO 28 I=1,MAXI
        WRITE(*,25) I
   25   FORMAT(//,1X,'Enter the number of coded data ',
     +               'values equal to ',I4,':  ')
        READ(*,*) NI(I)
   28 CONTINUE
C
      RETURN
      END
C
C
C----------------------------------------------------
C
       SUBROUTINE CONPLT(MAXI,NI,MUEST,SIGEST,LMAX)
C
       REAL*8      F(70,20),BD(9),XL(15),DUMF
       INTEGER     COUNT,PRNTFLG(74),FLAG(70,20),DFLAG
       CHARACTER*1 PRNTCHR(8)
      INTEGER      MAXI, NI(10)
      REAL*8       MUEST, SIGEST, LMAX,
     +             XMIN, XMAX, YMIN, YMAX, FMIN, FMAX,
     +             DX, DY, X, Y, X1, Y1,
     +             DF, P
C
        DATA PRNTCHR/'.','-','=','+','X','O','*',' '/
C
 29    WRITE(*,31)
 31    FORMAT(5(/),1X,'Enter the region of interest for the ',
     +                'contour plot',
     +/,1X,'(MU minimum, MU maximum, SIGMA minimum, SIGMA maximum)',
     +/,1X,'---->  ')
       READ(*,*) XMIN, XMAX, YMIN, YMAX
       IF (YMIN .LE. 0.0D0) YMIN=0.01D0
       WRITE(*,32)
 32    FORMAT(9(/))
C
       WRITE(*,*) ' Computing . . . please wait . . .'
       WRITE(*,32)
       FMIN=999999999.9D0
       FMAX=-999999999.9D0
       DX=(XMAX-XMIN)/64.0D0
       DY=(YMAX-YMIN)/16.0D0
C
       DO 10 J=1,17
        Y=YMAX-DBLE(J-1)*DY
        DO 20 I=1,65
         X=XMIN+DBLE(I-1)*DX
         CALL GETLOG(MAXI,NI,X,Y,DUMF,DFLAG)
         F(I,J) = DUMF
         FLAG(I,J) = DFLAG
         IF (FLAG(I,J) .EQ. 1) GOTO 20
         IF (F(I,J) .GT. FMAX) THEN
          FMAX=F(I,J)
          X1=X
          Y1=Y
         END IF
         IF (F(I,J) .GT. LMAX) THEN
          LMAX=F(I,J)
          MUEST=X
          SIGEST=Y
         END IF
         IF (F(I,J) .LT. FMIN) THEN
          FMIN=F(I,J)
C         X0=X
C         Y0=Y
         END IF
  20    CONTINUE
  10   CONTINUE
C
       DF=(FMAX-FMIN)/7.0D0
       BD(1)=FMIN
       BD(8)=FMAX
       DO 1 I=2,7
       BD(I)=FMIN + DF*DBLE(I-1)
  1    CONTINUE
       WRITE(*,2)
  2    FORMAT(///,' SIGMA',3X,'|',15X,'LOG L(MU, SIGMA)   ',
     +            ' . - = + X O * ')
       WRITE(*,3)
  3    FORMAT(9X,'|')
       COUNT=3
       DO 30 J=1,17
        DO 40 I=1,65
         IF (FLAG(I,J) .EQ. 0) GOTO 992
          PRNTFLG(I)=8
          GOTO 40
  992    P=F(I,J)
         IF (P .LT. BD(2)) THEN
          PRNTFLG(I)=1
         ELSE IF (P .LE. BD(3)) THEN
          PRNTFLG(I)=2
         ELSE IF (P .LE. BD(4)) THEN
          PRNTFLG(I)=3
         ELSE IF (P .LE. BD(5)) THEN
          PRNTFLG(I)=4
         ELSE IF (P .LE. BD(6)) THEN
          PRNTFLG(I)=5
         ELSE IF (P .LE. BD(7)) THEN
          PRNTFLG(I)=6
         ELSE
          PRNTFLG(I)=7
         END IF
 40     CONTINUE
C
        Y=YMAX-DY*DBLE(J-1)
        IF (COUNT .EQ. 3) THEN
         WRITE(*,4) Y,(PRNTCHR(PRNTFLG(K)),K=1,65)
  4      FORMAT(1X,F8.2,'+',65A1)
         COUNT=0
        ELSE
         WRITE(*,5) (PRNTCHR(PRNTFLG(K)),K=1,65)
  5      FORMAT(9X,'|',65A1)
         COUNT=COUNT+1
        END IF
 30    CONTINUE
C
C
C  Put the labels under the horizontal axis.
C
      WRITE (*,945)
  945 FORMAT (10X,'-+-------+-------+-------+----',
     +'---+-------+-------+-------+-------+----')
      DX=(XMAX-XMIN)/8.0D0
      DO 950 I=1,9
        XL(I)=XMIN+DBLE(I-1)*DX
  950 CONTINUE
C
      WRITE (*,960) (XL(I),I=1,9)
  960 FORMAT(7X,9(F7.3,1X),'MU')
C
      CALL PAUSE
C
       WRITE(*,15)
 15    FORMAT(//////,10X,'SYMBOL',15X,'RANGE OF F',/)
       DO 70 I=1,7
        WRITE(*,16) PRNTCHR(I),BD(I),BD(I+1)
 16     FORMAT(12X,A1,12X,F9.4,1X,'-',1X,F9.4)
 70    CONTINUE
       WRITE(*,13)
 13    FORMAT(2(/),10X,'A POINT OF MAXIMUM ON THIS PLOT IS:',/,7X,
     + 'MU          SIGMA     LN L(MU,SIGMA)',/)
       WRITE(*,100) X1,Y1,FMAX
C
       WRITE(*,12)
 12    FORMAT(/,10X,'THE HIGHEST POINT OF MAXIMUM FOUND '
     +              'IN THIS SESSION IS:',/,7X,
     + 'MU          SIGMA     LN L(MU,SIGMA)',/)
       WRITE(*,100) MUEST,SIGEST,LMAX
C
100    FORMAT(3X,F9.4,3X,F9.4,4X,F9.4)
       WRITE(*,14)
 14    FORMAT()
C
      CALL PAUSE
       RETURN
       END
C
C
C----------------------------------------------------
C
       SUBROUTINE BNDPLT(MAXI,NI,MUEST,SIGEST,LMAX)
C
       REAL*8      F(70,20),BD(9),CONF(7),XL(15),DUMF
       INTEGER     COUNT,PRNTFLG(74),FLAG(70,20),DFLAG
       CHARACTER*1 PRNTCHR(8)
      INTEGER      MAXI, NI(10)
      REAL*8       MUEST, SIGEST, LMAX,
     +             XMIN, XMAX, YMIN, YMAX,
     +             PT75, PT90, PT95, PT975, PT99, PT995, PT999,
     +             DX, DY, X, Y, P
C    +             DX, DY, X, Y, X1, Y1, X0, Y0,
C    +             DF, P
C
        DATA PRNTCHR/'.','-','=','+','X','O','*',' '/
        DATA CONF/99.9D0,99.5D0,99.0D0,97.5D0,95.0D0,90.0D0,75.0D0/
C
 29    WRITE(*,31)
 31    FORMAT(/,1X,'Enter the region of interest for the confidence ',
     +'band plot',
     +/,1X,'(MU Minimum, MU Maximum, SIGMA Minimum, SIGMA Maximum)',
     +/,1X,'---->  ')
       READ(*,*) XMIN, XMAX, YMIN, YMAX
       IF (YMIN .LE. 0.0D0) YMIN=0.01D0
       WRITE(*,32)
 32    FORMAT(9(/))
C
       WRITE(*,*) 'Computing . . . please wait . . .'
       WRITE(*,32)
       DX=(XMAX-XMIN)/64.0D0
       DY=(YMAX-YMIN)/16.0D0
C
      PT75  = LMAX - 0.50D0* 2.27D0
      PT90  = LMAX - 0.50D0* 4.61D0
      PT95  = LMAX - 0.50D0* 5.99D0
      PT975 = LMAX - 0.50D0* 7.38D0
      PT99  = LMAX - 0.50D0* 9.21D0
      PT995 = LMAX - 0.50D0*10.60D0
      PT999 = LMAX - 0.50D0*13.80D0
C
       DO 10 J=1,17
        Y=YMAX-DBLE(J-1)*DY
        DO 20 I=1,65
         X=XMIN+DBLE(I-1)*DX
         CALL GETLOG(MAXI,NI,X,Y,DUMF,DFLAG)
         F(I,J) = DUMF
         FLAG(I,J) = DFLAG
         IF (FLAG(I,J) .EQ. 1) GOTO 20
  20    CONTINUE
  10   CONTINUE
C
      BD(1)=PT999
      BD(2)=PT995
      BD(3)=PT99
      BD(4)=PT975
      BD(5)=PT95
      BD(6)=PT90
      BD(7)=PT75
      BD(8)=LMAX
C
       WRITE(*,2)
  2    FORMAT(///,' SIGMA',3X,'|',15X,'CONFIDENCE BANDS ',
     +            ' . - = + X O * ')
       WRITE(*,3) MUEST, SIGEST
  3    FORMAT(9X,'|',8X,'ESTIMATE OF MU = ',F8.3,
     +                  '     ESTIMATE OF SIGMA = ',F8.3)
       COUNT=3
       DO 30 J=1,17
        DO 40 I=1,65
         IF (FLAG(I,J) .EQ. 0) GOTO 992
          PRNTFLG(I)=8
          GOTO 40
  992    P=F(I,J)
         IF (P .LT. BD(1)) THEN
          PRNTFLG(I)=8
         ELSE IF (P .LE. BD(2)) THEN
          PRNTFLG(I)=1
         ELSE IF (P .LE. BD(3)) THEN
          PRNTFLG(I)=2
         ELSE IF (P .LE. BD(4)) THEN
          PRNTFLG(I)=3
         ELSE IF (P .LE. BD(5)) THEN
          PRNTFLG(I)=4
         ELSE IF (P .LE. BD(6)) THEN
          PRNTFLG(I)=5
         ELSE IF (P .LE. BD(7)) THEN
          PRNTFLG(I)=6
         ELSE
          PRNTFLG(I)=7
         END IF
 40     CONTINUE
C
        Y=YMAX-DY*DBLE(J-1)
        IF (COUNT .EQ. 3) THEN
         WRITE(*,4) Y,(PRNTCHR(PRNTFLG(K)),K=1,65)
  4      FORMAT(1X,F8.2,'+',65A1)
         COUNT=0
        ELSE
         WRITE(*,5) (PRNTCHR(PRNTFLG(K)),K=1,65)
  5      FORMAT(9X,'|',65A1)
         COUNT=COUNT+1
        END IF
 30    CONTINUE
C
C
C  Put the labels under the horizontal axis.
C
      WRITE (*,945)
  945 FORMAT (10X,'-+-------+-------+-------+----',
     +'---+-------+-------+-------+-------+----')
      DX=(XMAX-XMIN)/8.0D0
      DO 950 I=1,9
        XL(I)=XMIN+DBLE(I-1)*DX
  950 CONTINUE
C
      WRITE (*,960) (XL(I),I=1,9)
  960 FORMAT(7X,9(F7.3,1X),'MU')
C
      CALL PAUSE
C
       WRITE(*,15)
 15    FORMAT(12(/),10X,'SYMBOL',15X,'LEVEL OF CONFIDENCE',/)
       DO 70 I=1,7
        WRITE(*,16) PRNTCHR(I),CONF(I)
 16     FORMAT(12X,A1,26X,F5.2)
 70    CONTINUE
       WRITE(*,14)
 14    FORMAT(5(/))
C
      CALL PAUSE
       RETURN
       END
C


C--------------------------------------------------------------------------
C
       SUBROUTINE GETLOG(MAXI,NI,MU,SIGMA,LOGLIK,PFLAG)
C
       REAL*8  PHI(11), P(10), MU, SIGMA, LOGLIK
       INTEGER MAXI, NI(10), PFLAG
C
      CALL GETPHI(MAXI,MU,SIGMA,PHI)
C
      PFLAG = 0
      CALL GETP(MAXI,PHI,P,PFLAG)
      LOGLIK=0.0D0
      IF (PFLAG .EQ. 1) GOTO 120
C
      DO 100 I=1,MAXI
        LOGLIK = LOGLIK + NI(I) * DLOG ( P(I) )
  100 CONTINUE
C
  120 RETURN
      END
C
C
C--------------------------------------------------------------------
C
      SUBROUTINE GETPHI(MAXI,MU,SIGMA,PHI)
C
      REAL*8  MU, SIGMA, PHI(11), Z, NORMCP
      INTEGER MAXI, MAXP1
C
      MAXP1 = MAXI + 1
C
      DO 100 I=1,MAXP1
        Z = ((I - 0.5D0) - MU) / SIGMA
        PHI(I) = NORMCP(Z)
  100 CONTINUE
C
      RETURN
      END
C
C----------------------------------------------------
C
      DOUBLE PRECISION FUNCTION NORMCP(Z)
C
C  This function returns the standard normal cumulative
C  probability of being less than Z.
C
      REAL*8 Z,C(7),Y,Q,S
      DATA C/.319381530D0,-.356563782D0,1.781477937D0,
     +       -1.821255978D0,1.330274429D0,.2316419D0,
     +       2.506628725D0/
      IF (Z .GT. 6.0D0) THEN
        NORMCP = 1.0D0
C
      ELSEIF (Z .LT. -6.0D0) THEN
        NORMCP = 0.0D0
C
      ELSE
        Y=Z
        IF (Z .LT. 0.0D0) Y=-Z
        Q=1.0D0/(1.0D0+C(6)*Y)
        S=((((C(5)*Q+C(4))*Q+C(3))*Q+C(2))*Q+C(1))*Q
        NORMCP = S*DEXP(-Y*Y/2.0D0)/C(7)
        IF (Z .GT. 0.0D0) NORMCP=1.0D0-NORMCP
      ENDIF
C
  100 RETURN
      END
C
C-----------------------------------------------------
C
C
      SUBROUTINE GETP(MAXI,PHI,P,PFLAG)
C
      REAL*8  PHI(11), P(10)
      INTEGER MAXI, PFLAG
C
      P(1) = PHI(2) - PHI(1)
      IF (P(1) .LE. 0.0D0) PFLAG=1
C
      DO 10 I=2,MAXI
        P(I) = PHI(I+1) - PHI(I)
        IF (P(I) .LE. 0.0D0) PFLAG=1
   10 CONTINUE
C
      RETURN
      END
C
C
C----------------------------------------------------
C
      SUBROUTINE PAUSE
C
      CHARACTER*1 GO
C
C  Pause to let the user view a display.
C
      WRITE(*,70)
   70 FORMAT(/,1X,'STRIKE CARRIAGE RETURN ',
     +'WHEN YOU WISH TO CONTINUE  ')
      READ(*,75) GO
   75 FORMAT(A1)
C
      RETURN
      END
C
C----------------------------------------------------
C
      SUBROUTINE CHK(NUM)
C
      INTEGER NUM
C
      WRITE(*,*) 'CHECKPOINT ',NUM
C
      RETURN
      END
C
C----------------------------------------------------