Writing an Algorithm for maximum likelihood estimation in MATLAB

Question

I know that MATLAB has a package for maximum likelihood estimation, but for educational purposes I'm writing by myself an algorithm that gives me back estimates. Now, I have written a function that I'm trying to minimize (since I'm using the negative log likelihood). Here it is:

function ml = two_var(param, data)
mu = param(1);
sigma = param(2);
n=numel(data);
sumto = 0;
for i=1:n
    c = data(i)- mu;
    sumto = sumto + c;
    ml = n/2*log(2*pi)+n/2*log(sigma^2)+1/(2*sigma^2)*sumto^2;
end

This code concerns the estimation of a gaussian distibution. Now, the problem that I have is that this function does not seem to be a valid fminunc input... How can I circumvert the problem? What am I doing wrong? Thanks to anybody who wants to help ;)

why is it not a valid input? A way of circumventing the problem is writing your own minimization function instead of using matlabs inbuilt.... — Ander Biguri
Have you called fminunc with fminunc(@(param) two_var(param, data), param0 ) ? — Jommy
Yes... It tells me "failure in initial user-supplied objective function evaluation. FMINUNC cannot continue." — james42
Your loop does not really make sense. Take a look what is happening to sumto - on each step of the loop, you are adding c to it, and then computing again the variable ml (whose value at the final iteration is what your function will return). You should probably get rid of the loop altogether and just apply vectorized operations to data. Then use sum to compute the value of the likelihood function. — mikkola

james42 james42 · Accepted Answer · 2016-06-15T20:49:21

Thanks to everybody that helped in commenting the question. Well, I have used all your useful comments in order to improve my code. Here is my final result.

function ml = t_var(param)
rng default;
data = random('norm',0,1,[400,1]);
z = (data - param(1)) ./ param(2);
L = -.5.*z.*z - log(sqrt(2.*pi).*param(2));
ml = -sum(numel(data).*L);

This code is undoubtedly easier to read; moreover, it makes use of the fast vector operations that are allowed in MATLAB. After defining this function, is then easy to call it with

x0 =[0 0]
[x,fval] = fminsearch(@t_var,x0,optimset('TolX',1e-19,'MaxFunEvals',1000,'MaxIter',1000))

And obtain ML estimates that are consistent with the ones that you can get otherwise.

Thanks a lot to everybody who gave me some hints! ;)

Writing an Algorithm for maximum likelihood estimation in MATLAB

1 Answers