Maximizing likelihood not working

Hi, I am trying to maximize the likelihood of a sample using MATLAB. My issue is that, when testing my code against a result i know should return some values, it just returns the initial guess. So I changed many things but cannot get it to work properly, the same error happens again and again. I am a rather new user, so perhaps there is something easy I am missing. This is the code:

mb =  [60    20    20     0    60    20    20     0;...
       5    10    85     0     5    10    85     0;...
       0     5    60    35     0     5    60    35;...
       0     0     1    99     0     0     1    99];...

y=[3,1,1/3]'; %these are the values of xi/(1-xi) for the constraints.

pi=0.63; %this is the prior
pj=(1-pi)/pi;

%define likelihood function

lik=@(x)-(10^130)*(x(1,1)^(mb(1,1)+mb(1,5)))*(x(1,2)^(mb(1,2)+mb(1,6)))*(x(1,3)^(mb(1,3)+mb(1,7)))*(x(1,4)^(mb(1,4)+mb(1,8)))...
        *(x(2,1)^(mb(2,1)+mb(2,5)))*(x(2,2)^(mb(2,2)+mb(2,6)))*(x(2,3)^(mb(2,3)+mb(2,7)))*(x(2,4)^(mb(2,4)+mb(2,8)))...
        *(x(3,1)^(mb(3,1)+mb(3,5)))*(x(3,2)^(mb(3,2)+mb(3,6)))*(x(3,3)^(mb(3,3)+mb(3,7)))*(x(3,4)^(mb(3,4)+mb(3,8)))...
        *(x(4,1)^(mb(4,1)+mb(4,5)))*(x(4,2)^(mb(4,2)+mb(4,6)))*(x(4,3)^(mb(4,3)+mb(4,7)))*(x(4,4)^(mb(4,4)+mb(4,8)));

%initial values for iterations.

x0=[0.6 0.2 0.1 00.1;0.05 0.1 0.8 0.05;0.05 0.05 0.55 0.35;0.1 0.1 0.1 0.7];

%no inequality constraints
A=[];
b=[];
%make each line of x sum to 1
Aeq=[1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0; ...
     0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0; ...
     0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0; ...
     0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1];
beq=ones(4,1);

%lower bound zero, upper bound 1
lb=zeros(4,4);
ub=ones(4,4);

%non-linear constraints in cons_am
nonlcon=@cons_am;
options = optimoptions(@fmincon,'TolCon',1.0000e-09);
x=fmincon(lik,x0,A,b,Aeq,beq,lb,ub,nonlcon);

And the no linear constraints are:

function [c, ceq] = cons_am(x)
c(1)=3-0.5873*(0.03*x(1,1)+0.34*x(2,1)+0.46*x(3,1)+0.17*x(4,1))/(0*x(1,1)+0.14*x(2,1)+0.41*x(3,1)+0.45*x(4,1));
c(2)=  0.5873*(0.03*x(1,2)+0.34*x(2,2)+0.46*x(3,2)+0.17*x(4,2))/(0*x(1,2)+0.14*x(2,2)+0.41*x(3,2)+0.45*x(4,2))-3;
c(3)=1-0.5873*(0.03*x(1,2)+0.34*x(2,2)+0.46*x(3,2)+0.17*x(4,2))/(0*x(1,2)+0.14*x(2,2)+0.41*x(3,2)+0.45*x(4,2));
c(4)=  0.5873*(0.03*x(1,3)+0.34*x(2,3)+0.46*x(3,3)+0.17*x(4,3))/(0*x(1,3)+0.14*x(2,3)+0.41*x(3,3)+0.45*x(4,3))-1;
c(5)=(1/3)-0.5873*(0.03*x(1,3)+0.34*x(2,3)+0.46*x(3,3)+0.17*x(4,3))/(0*x(1,3)+0.14*x(2,3)+0.41*x(3,3)+0.45*x(4,3));
c(6)=  0.5873*(0.03*x(1,4)+0.34*x(2,4)+0.46*x(3,4)+0.17*x(4,4))/(0*x(1,4)+0.14*x(2,4)+0.41*x(3,4)+0.45*x(4,4))-(1/3);
ceq=[];

The non linear restrictions have to do with the probabilities being assigned to some intervals, the others have to do with adding up to one, etc. Using the 'mb' matrix as is defined should return the same values but each divided by 10 (did the calculations in Excel). The (10^130) in the beginning of the function is to avoid the numbers becoming super-small.

Thanks very much in advance.

Answer 1

[Alan Weiss]

I think that you would have more success by maximizing the log of the likelihood, which should avoid many numerical issues. And if you want fmincon to use you options, you need to pass the options in:

x = fmincon(lik,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);

Alan Weiss

MATLAB mathematical toolbox documentation

[Juan Gambetta]

The log is a very nice idea. Also I corrected the options argument. Thank you very much, will check if it works.