Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: fmincon stuck...
Replies: 1   Last Post: Jul 17, 2009 7:19 PM

 Messages: [ Previous | Next ]
 Andrea Asoni Posts: 10 Registered: 4/15/09
fmincon stuck...
Posted: Jul 17, 2009 5:10 PM

Hi,

I am a matlab beginner... I hope this question is not too dumb. I would like to solve a constrained maximization problem. I am using the fmincon function on Matlab. I have the following problem:

fmincon does not move away from initial condition. every iteration of the optimization process gives me the same value for the function and the parameters. Any suggestion? Is there anything I should be thinking of except an error in my code?

% MAIN PROGRAM

data=[stat2 stat23 stat24 stat3 stat32 stat34 stat4 stat42 stat43 id co age nam sex marr kid urb locus ld35 ldsq35];

% INITIAL CONDITIONS
tes0=[0.2 1]; ten0=[0.2 1]; tse0=[0.2 1]; tsn0=[0.2 1]; tne0=[0.2 1]; tns0=[0.2 1];
bes0=[0.02 0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 0.02 -0.02]; ben0=[0.02 0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 0.02 -0.02];
bse0=[0.02 0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 0.02 -0.02]; bsn0=[0.02 0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 0.02 -0.02];
bne0=[0.02 0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 0.02 -0.02]; bns0=[0.02 0.02 0.02 -0.02 0.02 -0.02 0.02 -0.02 0.02 0.02 -0.02];

t0=[tes0 ten0 tse0 tsn0 tne0 tns0];
b0=[bes0 ben0 bse0 bsn0 bne0 bns0];

% MATRIX WITH ALL COMBINATION OF THETAS
M=allcomb(tes0, ten0, tse0, tsn0, tne0, tns0);
M=M';
sla=size(M);
sla=sla(2);

% PROBABILITES
p0=1/sla*ones(1,sla);

% INITIAL CONDITIONS AND PARAMETERS
init=[t0,p0,b0];
t=t0;
b=b0;
p=p0;
z=[t p b];

% HANDLES
fhandle=@(z) jmplik(z,data);
mycons=@(z) jmpcons(z,data);

% OPTIMIZATION
[z,fval]=fmincon(fhandle,init,[],[],[],[],[],[],mycons)

----------------- MAIN FUNCTION
function L=jmplik(z,data);

tes=z(1,(1:2)); ten=z(1,(3:4)); tse=z(1,(5:6)); tsn=z(1,(7:8)); tne=z(1,(9:10)); tns=z(1,(11:12)); p=z(1,(13:76));
bes=z(1,(77:87)); ben=z(1,(88:98)); bse=z(1,(99:109)); bsn=z(1,(110:120)); bne=z(1,(121:131)); bns=z(1,(132:142));

stat2=data(:,1); stat23=data(:,2); stat24=data(:,3); stat3=data(:,4); stat32=data(:,5); stat34=data(:,6); stat4=data(:,7); stat42=data(:,8); stat43=data(:,9);
id=data(:,10); co=data(:,11); age=data(:,12); nam=data(:,13); sex=data(:,14); marr=data(:,15); kid=data(:,16); urb=data(:,17); locus=data(:,18); ld35=data(:,19); ldsq35=data(:,20);

n=length(id);
const=ones(n,1);

% MATRIX WITH ALL COMBINATION OF THETAS
M=allcomb(tes, ten, tse, tsn, tne, tns);
M=M';
sla=size(M);
sla=sla(2);

% DEFINE LAMBDAS
for j=1:sla;
yes(:,j)=bes(1).*const + bes(2).*co + bes(3).*age + bes(4).*nam + bes(5).*sex + bes(6).*marr + bes(7).*kid + bes(8).*urb + bes(9).*locus + bes(10).*ld35 + bes(11).*ldsq35 + M(1,j);
yen(:,j)=similar formula...
yse(:,j)=similar formula...
ysn(:,j)=similar formula...
yne(:,j)=sim for
yns(:,j)=sim for
end

laes=1.+ exp(-yes);
laen=1.+ exp(-yen);
lase=1.+ exp(-yse);
lasn=1.+ exp(-ysn);
lane=1.+ exp(-yne);
lans=1.+ exp(-yns);

laes=1./laes;
laen=1./laen;
lase=1./lase;
lasn=1./lasn;
lane=1./lane;
lans=1./lans;

% CONSTRUCT HAZARD RATES AS MULTIPLICATION OF LAMBDAS
for j=1:sla;
hr(:,j)=stat2.*((1.-laes(:,j)).*(1.-laen(:,j))) + stat23.*(laes(:,j))+ stat24.*(laen(:,j))+...
stat3.*((1.-lase(:,j)).*(1.-lasn(:,j))) + stat32.*(lase(:,j))+ stat34.*(lasn(:,j))+...
stat4.*((1.-lane(:,j)).*(1.-lans(:,j))) + stat42.*(lane(:,j))+ stat43.*(lans(:,j));
end

% MULTIPLY THE HAZARD RATES FOR EACH ID: CALCULATE f(t|theta)
test = [id hr]; nrcol=size(test); nrcol=nrcol(2);
ID=unique(test(:,1)); len=length(ID);
prov=zeros(len,nrcol-1);

for i=1:len
gg=ID(i);
ind=find(test(:,1)==gg);
temp=test(ind,(2:nrcol));
prov(i,:)=prod(temp);
end;

% WEIGHTED SUM OF f(t|thetas) OVER THE POSSIBLE REALIZATIONS OF THETAS
p=p';
fsum=prov*p;

% TAKE LOGS AND SUM OVER THE LOGS TO GET THE LIKELIHOOD
fsum=log(fsum);
L=sum(fsum);
L=-L
z
pause

---FUNCTION FOR COSTRAINTS
function [c, ceq]=jmpcons(z,data)

p=z(1,(13:76));
sp=length(p);

c=-p;

ceq=p*ones(sp,1)-1;

Date Subject Author
7/17/09 Andrea Asoni
7/17/09 Matt