I am trying to teach myself maximum likelihood estimation using the Newton-Raphson method and related iterative methods. I don't understand the link between the hessian, the expected value of the hessian, the information matrix, and the expected information matrix.

Most resources say that you can compute estimates of parameters using the inverse of the hessian. (I have been referring to this pdf on maximum likelihood.)

I am testing the method by generating normally distributed data and trying to estimate the mean and variance. When I use the hessian method the variance seem to diverge. If I use a closed-form solution for the expected information matrix, the iterative method works. I want to use this method for more complex expressions. I figured it would be easier to stick with the hessian. How do I compute the expected value of the hessian?

Here is my Matlab code:

clear all; close all;
N = 500;
mu = 1; var = 1;
% Generate normally distributed data with mean mu and variance var
data = mu + sqrt(var)*randn(N,1);

Hess = zeros(2);
I = zeros(2);
estim = [1.4; 1.4];

% Closed-form estimates
MLE_mean = mean(data)
MLE_var = sum((data-MLE_mean).^2)/N

for iter = 1:400
 m = estim(1);
 v = estim(2);

 dldm = (sum(data) - N*m)/v;             % partial derivative wrt mean

 dldv = (sum((data-m).^2)/v  - N)/(2*v); % partial derivative wrt variance

 grad = [dldm; dldv];                    %gradient of log-likelihood

 Hess(1,1) = -N/v;
 Hess(2,2) = -(sum(data) - N*m)/(v.^3) + N/(2*v.^2);

 Hess(1,2) = -(sum(data) - N*m)/(v.^2);
 Hess(2,1) = Hess(1,2);

 Inverse_of_Expected_Info(1,1) = v/N;
 Inverse_of_Expected_Info(2,2) = 2*(v^2)/N;

 %estim = estim + Inverse_of_Expected_Info*grad;   % This method works

 estim = estim -(pinv((Hess))*grad);            % This method doesn't work


Iterative_estimate_mean = estim(1)
Iterative_estimate_var = estim(2)

You are using unit step-length for your optimization, that's why your method diverges after 40-50 iterations. Here is the plot of your intermediate solutions: enter image description here

You should perform line search (e.g. backtracking) along your Newton directions.

  • $\begingroup$ Thank you! So the general method is okay then? $\endgroup$ Sep 23 '14 at 20:09
  • $\begingroup$ yourwelcome, yep $\endgroup$
    – Arya
    Sep 23 '14 at 21:49

I think, $Hess(2,2) = -\sum(data.*data) -2*m*\sum(data)+ \frac{N*m*m}{v.^3} + \frac{N}{2*v.^2}$

as per the MLE method. On substituting this value, I'm able to get the correct solution even with the pinv() command.


This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .