function [mle, info, n] = gammanr(x, theta0, TOL, Nmax)
% *************************************************************************
%
% FILE: gammanr.m
%
% Date Created: Feb 20, 2004
% Author: Mark E. Irwin
%
% Contents: Finds MLE and information for gamma distribution by Newton-Raphson
%
% Revision History
% Date          Name            Changes/Reasons
%
% *************************************************************************
%
% Input:
%
%   x: gamma observation vector
%   theta0: initial guess for MLE
%   TOL: convergence tolerance
%   Nmax: number of functional iterations
%
% Output:
%
%   mle: mle
%   info: information,
%   n: number of iterations until convergence
%
% *************************************************************************

% Initialize loop

score = zeros(2,1);
hess = zeros(2,2);

thetai = theta0;
i = 0;
diff = 2 * TOL;  % guarentees at least one pass through loop

n = length(x);
sumlogx = sum(log(x));
sumx = sum(x);

while (i <= Nmax && diff > TOL)
  i = i + 1;
  thetaold = thetai;
  lambda = thetaold(1);
  k = thetaold(2);
  
  % calculate score vector
  score(1) = (sumx / lambda - n * k) / lambda;
  score(2) = sumlogx - n * (log(lambda) + psi(k));
  
  % calculate hessian = - observed information
  hess(1,1) = (n * k - 2 * sumx / lambda) / lambda^2;
  hess(1,2) = -n / lambda;
  hess(2,1) = hess(1,2);
  hess(2,2) = -n * psi(1,k);
  
  thetai = thetaold - hess \ score;
  diff = max(abs(thetaold - thetai));
end

if (diff > TOL)
  warning('Method did not converage after Nmax steps');
end

% update hessian for final information matrix
lambda = thetai(1);
k = thetai(2);
hess(1,1) = (n * k - 2 * sumx / lambda) / lambda^2 ;
hess(1,2) = -n / lambda;
hess(2,1) = hess(1,2);
hess(2,2) = -n * psi(1,k);


mle = thetai;
info = -hess;
n = i;