# Problem with fmincon and many starting values

12 views (last 30 days)
MRC on 13 Feb 2014
Edited: Matt J on 13 Feb 2014
Hi all, I'm trying to code the following minimization problem using fmincon and many starting values but I get error messages because the objective function has the log function inside. How can I solve this problem? Here the code. Thank you!
n=2;
thetatrue=[1 1 2 2 0.2];
mu = [0 0];
sd = [1 thetatrue(2*n+1); thetatrue(2*n+1) 1];
X=data(:,1:n);
Y=data(:,n+1:size(data,2));
A1(:,2)=-X(:,1);
A1(:,1)=-1;
A2(:,2)=-X(:,2);
A2(:,1)=-1;
W1=(all(bsxfun(@eq,Y,[0 0]),2));
W2=1-W1;
cdfun=@(x) mvncdf( [A1*[x(1);x(3)], A2*[x(2);x(4)] ]...
,mu,[1 x(5); x(5) 1]);
options=optimset('Algorithm','interior-...
point','Display','iter','MaxIter',10000,'TolX',10^-30,'TolFun',10^-30);
% Grid of starting values for fmincon
alpha1=normrnd(thetatrue(1),1,1,5); %random draws from N(alpha1true,1)
alpha2=normrnd(thetatrue(2),1,1,5); %random draws from N(alpha2true,1)
beta1=normrnd(thetatrue(3),1,1,5); %random draws from N(beta1true,1)
beta2=normrnd(thetatrue(4),1,1,5); %random draws from N(beta2true,1)
rho=(-0.9:0.1:0.9);
CC = {alpha1,alpha2,beta1,beta2,rho};
b = cellfun(@numel,CC);
n1 = cumsum(b);
a = fullfact(b);
idx = bsxfun(@plus,a,[0, n1(1:end-1)]);
v = cat(2,CC{:});
theta0grid = v(idx); %grid of starting values
L_nstar=zeros(size(theta0grid,1),1);
Exitflag=zeros(size(theta0grid,1),1);
for j=1:size(theta0grid,1) %for each starting value
theta0=theta0grid(j,:);
[theta,fval,exitflag,output]=...
fmincon(@(x) log_lik(x,cdfun,W1,W2),theta0,[],[],[],[],[-Inf; -Inf; Inf; -Inf;-0.999],[+Inf; +Inf; +Inf; +Inf; 0.999],[],options);
L_nstar(j)=-fval;
Exitflag(j)=exitflag;
end
L_nstar=max(L_nstar);
and
function val=log_lik(theta,cdfun,W1,W2)
z=cdfun(theta);
val=-sum(W1.*log(z)+W2.*log(1-z));
MRC on 13 Feb 2014
for example when j=1, I get the following error
Error using barrier (line 22) Objective function is undefined at initial point. Fmincon cannot continue.
Error in fmincon (line 905) [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = barrier(funfcn,X,A,B,Aeq,Beq,l,u,confcn,options.HessFcn, ...

Matt J on 13 Feb 2014
Edited: Matt J on 13 Feb 2014
The sqp algorithm can recover from NaN's and Infs returned by the objective function and constraints.
Matt J on 13 Feb 2014
I think you need a starting point that is feasible, so that they have a foothold in the feasible set.
You should probably take a closer look at your initial values and at the values of mu and sigma to see why you're getting z=0 or 1. Since z=cdfun(theta) is from a Gaussian CDF, z should always assume values strictly between 0 and 1 where log(z) and log(1-z) should not be producing undefined values. Maybe you have overflow, in which case you're choosing highly extreme theta values, or you have sigma=0 somewhere...