!     Last change:  DOS  30 Jul 2000    1:59 pm
!                                       *** copyright 2000 ***
! *** filename chex97.f95               *** John F. Monahan **
!                                       **********************
      program chex97
!  Example 9.7 -- Cox's logistic regression problem
!   y(i) = number of ingots not ready for rolling
!   x(i) = heating times
!   m(i) = number of ingots with x(i) heating time
!
      implicit none
      real, dimension(4) :: x
      real b0,b1,w,rll,del0,del1,h00,h01,h11,s0,s1,det,gam,p,omp
      real flogit
      integer i,j
      integer, parameter :: n = 4
      integer, dimension(4) :: y,m
!                          put the data here
      data x/ 7., 14., 27., 51. /
      data m/ 55, 157, 159, 16 /
      data y/ 0, 2, 7, 3 /
!                          formats
  21  format(/'Iteration',i4,'   Betas',2f9.4,'   Log-like',f12.6)
  22  format(' Gradient Vector',3f12.6)
  23  format(' Hessian Matrix',3f15.4)
  24  format(' Newton/Scoring Step ',2f12.6)
  25  format(//' Covariance Matrix '/2x,f12.6/2x,2f12.6)
!
      OPEN( UNIT=6, FILE='chex97.out' )
!                          initial values -- no effect
      b0 = log(12./365.)
      b1 = 0.
!                          iterate using irwls
      do j = 1,10             ! up to 10 iterations
!                          initialize
      rll = 0.
      del0 = 0.
      del1 = 0.
      h00 = 0.
      h01 = 0.
      h11 = 0.
!                          go through the data n = 4
      do i = 1,n
      gam = b0 + b1*x(i)
!                          get probabilities
      p = flogit(gam)
      omp = flogit(-gam)
!                          log-likelihood
      rll = rll + y(i)*gam + m(i)*log(omp)
!                          gradient
      del0 = del0 + y(i) - m(i)*p
      del1 = del1 + (y(i)-m(i)*p)*x(i)
!                          hessian matrix
      w = m(i)*p*omp
      h00 = h00 - w
      h01 = h01 - w*x(i)
      h11 = h11 - w*x(i)*x(i)
      end do  ! loop on i
!                          write out some things
      write(6,21) j,b0,b1,rll
      write(6,22) del0,del1
      write(6,23) h00,h01,h11
!                          compute step
      det = h00*h11 - h01*h01
      s0 = (h11*del0 - h01*del1) / det
      s1 = (h00*del1 - h01*del0) / det 
      write(6,24) s0,s1
!                          move
      b0 = b0 - s0
      b1 = b1 - s1 
!                          end of iteration
      end do  ! loop on j
!                          compute covariance matrix
!                            watch signs and order
      h00 = - h00/det
      h11 = - h11 / det
      h01 =  h01 / det 
      write(6,25) h11,h01,h00
      stop
      end program chex97
      REAL FUNCTION FLOGIT(S)
!  COMPUTE THE DISTRIBUTION FUNCTION OF LOGISTIC
      IMPLICIT NONE
      REAL, INTENT(IN) :: S
      REAL AS,E
      AS = ABS(S)                                                      
      E = 0.D0                                                          
      IF(AS .LT. 125. ) E = EXP(-AS)                                 
      IF( S .LT. 0. ) FLOGIT = E/(1. + E)                           
      IF( S .GE. 0. ) FLOGIT = 1./(1. + E)                        
      RETURN                                                            
      END FUNCTION FLOGIT
! *** end of filename chex97.f95        *********************