Maximum Likelihood - Multinomial Probit Model

Dear all,

I am trying to implement an estimation of a multinomial probit model by maximization of the log likelihood function implemented in the following code:

    f=0;

    n=size(y,1);
    X=zeros(n,1);
    for i=2:size(y,1)

       for j=1:size(x,2)
           X(i)=X(i)+b(j+4)*x(i,j);
       end

       if y(i)==1
           P(i)=normcdf((b(1))+X(i)),0,1);
       elseif y(i)==2
           P(i)=normcdf((b(2)+X(i)),0,1)-normcdf((b(1)+X(i)),0,1);
       elseif y(i)==3
           P(i)=normcdf((b(3)+X(i)),0,1)-normcdf((b(2)+X(i)),0,1);
       elseif y(i)==4
           P(i)=normcdf((b(4)+X(i)),0,1)-normcdf((b(3)+X(i)),0,1);
       else
           P(i)=1-normcdf((b(4)+X(i)),0,1);
       end

       if P(i)~=0
         f=f-log(P(i));
       end

Although I know it is possible to fit this model using the mnrfit function (from statistics toolbox). I want to implement that because my I want to be able to modify the log likelihood function to implement custom functions.

For some reason I am not able to reproduce the results generated by mnrfit function. I am using fminsearch to perform the maximization.

Is there anything wrong?

Sincerely yours, Renoir

Answer 1

[Tom Lane]

Here's a simplification of your code along with what I think is the corresponding code to fit using mnrfit. I think if you tightened up the convergence limits, fminsearch would eventually produce the same answers as mnrfit, but I didn't push that too far.

I'd avoid conditionally evaluating log(P). If you wind up with P=0 (or worse, P<0), you want fminsearch to get a bad result so it can back away and try something else. You don't want to skip over the bad contributions.

function doit
rng(1)
x = randn(1000,3);
b = [-1 -.5 .5 1   -1 1 2]';

% generate random y using this x and b
xb = x*b(5:7);
xb = bsxfun(@plus,b(1:4)',xb);
p = normcdf(xb);
p(:,2:4) = diff(p,1,2);
p(:,5) = 1-sum(p,2);
y = mnrnd(1,p);

% fit using mnrfit
b1 = mnrfit(x,y,'model','ordinal','link','probit')

% convert to y category numbers, then fit using fminsearch  
ynum = y*(1:5)';
b2 = fminsearch(@(b) nlogf(b,ynum,x), [.1 .2 .3 .4 0 0 0])

function f=nlogf(b,y,x)
f=0;
    for i=1:size(y,1)
       X=0;
       for j=1:size(x,2)
           X=X+b(j+4)*x(i,j);
       end
       if y(i)==1
           P=normcdf(b(1)+X);
       elseif y(i)==2
           P=normcdf(b(2)+X)-normcdf(b(1)+X);
       elseif y(i)==3
           P=normcdf(b(3)+X)-normcdf(b(2)+X);
       elseif y(i)==4
           P=normcdf(b(4)+X)-normcdf(b(3)+X);
       else
           P=1-normcdf(b(4)+X);
       end
       f=f-log(max(0,P));
    end