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).