"Dev Vencappa" <firstname.lastname@example.org> wrote in message <email@example.com>... > Dear Matlab Users, > I am new user to Matlab, so apologies in advance if my question below appears too basic to regular users. > > I have to eventually write and maximise a customised log likelihood function in Matlab and as a starting point, decided to learn how to do this in Matlab by coding the log-likelihood functions for the exponential model with covariates (and right censoring). To begin with, my function script is coded as below and saved in a file called expcovariate_mle.m (Admittedly, the most advanced users of Matlab will recognise that the codes below may not be the most efficient way of solving this problem, and comments are most welcome. I could not manage code this in fminunc by the way, as I was repeatedly getting messages I could not relate to my codes, so any suggestion of doing the below using fminunc if possible would be very much appreciated) > > function like=expcovariate_mle(params,t,c,xb); > theta=exp(xb*params'); > like = -sum(c.*log(theta)-t.*theta); > grad= -sum(c./theta-t); > hess= -sum(-c./(theta.^2)); > end > > In my main script, I then code and evaluate the following: > > clear all; > warning off all > format long > > cd 'E:\research\' > xlsread('final_data.xls','Sheet1','A2:Y759'); > c=ans(:,8); > t=ans(:,1); > > data=ans(:,:); > settled=data(:,8); > def=data(:,17); > est=data(:,18); > lia=data(:,19); > nol=data(:,20); > pat=data(:,21); > lap=data(:,22); > nla=data(:,23); > > xb=[ones(758,1),def,est,lia,nol,pat,lap,nla]; > > options = optimset('MaxFunEvals',5000); > params0=ones(1,8); > [beta]=fminsearch(@(params) expcovariate_mle(params,t,c,xb),params0,options); > > It all runs fine, and returns values of beta, except that these are completely different from what I obtain using other software (STATA/LIMDEP). I appreciate that different software use different algorithms for this maximization problem. So I decided to test this incrementally by starting with a limited number of variables in the matrix xb, and comparing this with the results from the other software using the limited number of variables. So for example, in the above, I would initially start off with say: > > xb=[ones(758,1),def,est]; > params0=ones(1,3); > > The beta that matlab returns are exactly the same as what STATA/Limdep returns. Then I tried: > > xb=[ones(758,1),def,est,lia]; > params0=ones(1,4); > > and get the same results as STATA/Limdep. > > I get exactly the same results up to using 5 variables out of a list of 7: > > xb=[ones(758,1),def,est,lia,nol,pat]; > params0=ones(1,6); > > However, as soon as I introduce lap (or lap and nla) into the xb matrix, I get completely different results from what STATA/Limdep produces with the extra variables, i.e. the below gives completely different results from what STATA produces. > > xb=[ones(758,1),def,est,lia,nol,pat,lap]; > params0=ones(1,7); > > So I am not sure where I am going wrong with this. Any help would be very much appreciated. > > Regards, > Dev
How different are thes results? Did you try changing the tolerances using optimset?
I have my own nelder mead algorithm in Matlab that is based off of numerical methods in fortran 95 and always get slightly different results compared to fminsearch. Both tend to be acceptable.
Did you plot predicted vs actual and if so does this look reasonable?
I don't think it is surprising that you get different answers as you add more parameters. Try the tolerances. I have seen 1e-20 as extreme cases along with different initial guesses.
Without seeing your data, it's hard to know.
As a new Matlab user, you have definitely made nice progress.