Maximum Likelihood Estimation

In this post, i want to estimate the maximum likelihood by numerical solution using matlab / octave. numerical solution used is based on Newton’s method and central difference for evaluation of the derivative values of the loglikelihood.

here i will attach the matlab code which is simple and without built-in functions.

Assume we have 2 normally distributed independent variable.

let’s simulate them:

rng(123);                    % for reproducing the same simulated data
x1 = 2 + 5*randn(100,1);     % x1 is with mean = 2 and std = 5
x2 = 4 + 2*randn(100,1);     % x2 is with mean = 4 and std = 2
x = [ones(100,1) , x1 , x2]; % x is 100 * 3 matrix
y = 3 - 4.2 * x1 + 7.5 * x2 + randn(100,1); % dependent

In the above equation the dependent variable includes normally distributed random noise with mean 0 and standard deviation = 1.

Using the simulated data we can do linear regression to estimate the coefficients.

The negative log likelihood for normal distribution is as below:

function L = nloglikely(Y,X,b,sigma)
n=length(Y);
L = -(-n/2*log(2*pi*(sigma).^2) - ...
    1/2 * (Y - X*b)' * (Y-X*b) ./((sigma)^2));

how to calculate the maximum likelihood using above function. One way is to take the derivative of the log likelihood in respect to b (coefficients) and make it equal zero and estimate the coefficients solving the system of equation or ordinary least square method could be used to estimate the coefficients from the design matrix x, using b = inv(x’*x)*x’*y.

Anyway this solution work in linear regression, what if we want to model a generalized linear model with a exponential family distribution.

Here in this post i suggest a simple and numerical method for computation of maximum likelihood using:

function [b,dP,loglikelihoodE] = finitedifference(b,Y,X,MaxIter,Tol)
n = length(Y);
h = 0.01;
p = length(b);
b0 = zeros(size(b));
while sum(abs(b - b0)) >= Tol
c = zeros(size(b));b0 = b;
for o = 1:length(b)
 b1 = b;b2 = b;b_2=b;b_1=b;
 for k = 1:MaxIter
 b2(o) = b(o) + 2*h;
 b_2(o) = b(o) - 2*h;
 b1(o) = b(o) + h;
 b_1(o) = b(o) - h;
 dP = 0;
 for i = 1:length(Y)
 dP = dP + (Y(i) - X(i,:)*b).^2 ./ (n - p);
 end
 f(o) = (-L(Y,X,b2,dP,n) + 8*L(Y,X,b1,dP,n) - 8*L(Y,X,b_1,dP,n) + L(Y,X,b_2,dP,n))/12/h; 
 ff(o) = (-L(Y,X,b2,dP,n) + 16*L(Y,X,b1,dP,n) - 30*L(Y,X,b,dP,n) + 16*L(Y,X,b_1,dP,n) - L(Y,X,b_2,dP,n))/12/h^2; 
 c(o) = b(o) - f(o)/ff(o);
 if abs(b(o) - c(o)) < Tol
 break;
 end
 b(o) = c(o);
 end
end
dP
b = c
loglikelihoodE = -L(Y,X,b,dP,n)
end

function l=L(Y,X,b,sigma,n)
l = -(-n/2*log(2*pi*(sigma).^2) - 1/2 * (Y - X*b)' * (Y-X*b) ./((sigma)^2));
end
end

This function use central difference O(4) to calculate the derivative value of the negative log likelihood and second derivative value of negative log likelihood function. Then to calculate the maximum likelihood or minimum of the negative loglikelihood it will use Newton’s method. for this simulated data running the above function returns:

> finitedifference([27,22,57]',y,x,10000,.00001)

dP =

1.3676

b =

2.5374
 -4.1897
 7.5896

loglikelihoodE =

-158.6637

now let’s compare this results with the built-in matlab function fitglm.

> glm = fitglm(x,y, 'intercept', false)
Generalized Linear regression model:
 y ~ 1 + x1 + x2 + x3
 Distribution = Normal

Estimated Coefficients:
               Estimate     SE    tStat     pValue 
               ________ ________ _______ ___________

(Intercept)    2.5374     0.2778  9.1338   1.08e-14
 x2           -4.1897   0.029621 -141.45   3.1832e-113
 x3            7.5896   0.062166  122.09   4.0387e-107

100 observations, 97 error degrees of freedom
Estimated Dispersion: 1.37
F-statistic vs. constant model: 1.63e+04, p-value = 2.48e-123

> glm.LogLikelihood

-156.0242

As we see in the above results, The estimated dispersion parameter 1.37 is same in both models and coefficients are same. the maximum likelihood is also close using both methods. (-158 and -156).

 

 

Leave a Reply

Your email address will not be published. Required fields are marked *