# chex98.s                         J F Monahan, June 2001
 # IRWLS for Finney logistic regression example
 # Model has y as binomial( m, p ) with
 nobs <- 39
 p <- 3
 X <- scan("finney.dat")
 X <- matrix( X, p, nobs)
 y <- X[3,]
 y
 X <- matrix( c(rep(1,nobs),log(X[1,]),log(X[2,])),nobs, 3 )
 X 
 # initial values
 b <- c( log(sum(y)/(nobs-sum(y))), 0 , 0)
 m <- rep(1,nobs)
 for (j in 1:10) {          # ten iterations
 cat("***************** iteration",j,"\n")
 print('beta')
 print(b)
 gam <- X %*% b
 p <- exp(gam)/(1+exp(gam))
 omp <- 1/(1+exp(gam))
 rll <- sum(y*gam) + sum(m*log(omp))
 print("log-like")
 print(rll)
 del <- t(X) %*% (y-m*p)
 print("gradient")
 print(del)
 w <- m*p*omp
 h <- - t(X) %*% diag(w[,1],nrow=length(m)) %*% X
 print("Hessian")
 print(h)
 s <- solve(h,del)
 print("Newton/Scoring Step")
 print(s) 
 b <- b - s                 # update
                 }          # end iteration loop on j
 print("Covariance Matrix")
 solve(-h)
 rm(list=ls())