/* Interatively Reweighted Least Squares with Poisson Errors

Author: Ted Thompson   tat5@ccdddt1.em.cdc.gov
        Centers for Disease Control and Prevention
Date: 18 Jul 95
Provided without guarantee for public non-commercial use.


The following procedure fits a generalized linear model via iteratively
reweighted LS with log link and poisson errors.
The inputs are y , the response variable, x, the predictors,
and off, an offset.   The outputs are i, number of iterations,
dev, deviance value, beta, estimmated coefficients, and var,
the covariance matrix of the beta values.
*/

proc( 4 ) = poi( y, x, off );
local beta, beta0, c, n, p, z, eps, i, var, dev, xp;

beta0 = zeros( cols(x), 1 );
beta = ones( cols(x), 1 );
eps = 1.0e-8;
i = 0;

do while (abs(beta-beta0) > eps) and (i <= 50);
  beta0 = beta;
    if i eq 0;
      n = ln(y+0.5);
    else;
      n = x * beta0 + off;
    endif;
  p = exp(n);
  z = n + ( (y-p) ./ p ) - off;
  xp = x .* p;
  c = chol( xp'x );
  beta = cholsol( xp'z, c );
i = i + 1;
endo;

var = cholsol( eye(rows(c)), c );
dev = 2 * sumc( y.*ln(y./p) - (y-p) );

retp( i, dev, beta, var );
endp;

