// Newton-Raphson for Poisson log-likelihood
clear all
use accident3

mata:
real scalar pll(real vector y, real matrix X, real vector b)
{
    real vector  xb
    xb = X*b
    return(sum(-exp(xb) + y:*xb - lnfactorial(y)))
}

void GetDerives(real vector y, real matrix X, real vector theta, g, Hi)
{
	real vector exb

	exb = exp(X*theta)
	g   =  (quadcolsum((y - exb):*X))'
	Hi  = quadcross(X, exb, X)
	Hi = -1*cholinv(Hi)
}

real vector tupdate(                 ///
	real scalar lambda,          ///
	real vector theta_s,          ///
	real vector g_s,              ///
	real matrix Hi_s)
{
	return (theta_s - lambda*Hi_s*g_s)
}

real vector GetUpdate(            ///
    real vector y,                ///
    real matrix X,                ///
    real vector theta_s,          ///
    real vector g_s,              ///
    real matrix Hi_s)
{
    real scalar lambda 
    real vector theta_s1

    lambda = 1
    theta_s1 = tupdate(lambda, theta_s, g_s, Hi_s)
    while ( pll(y, X, theta_s1) <= pll(y, X, theta_s) ) {
        lambda   = lambda/2
        if (lambda < 1e-11) {
            printf("{red}Cannot find parameters that produce an increase.\n")
            exit(error(3360))
        }
        theta_s1 = tupdate(lambda, theta_s, g_s, Hi_s)
    }
    return(theta_s1)
}


y = st_data(., "accidents")
X = st_data(., "cvalue kids traffic")
X = X,J(rows(X), 1, 1)

b  =  J(cols(X), 1, .01)
GetDerives(y, X, b, g=., Hi=.)
gz = .
while (abs(gz) > 1e-11) {
	bs1 = GetUpdate(y, X, b, g, Hi)
	b   = bs1
	GetDerives(y, X, b, g, Hi)
	gz = g'*Hi*g
	printf("gz is now %8.7g\n", gz)
}
printf("Converged value of beta is\n")
b

end