Minimize function by choosing a variable that enters indirectly

Hello,

I am running a simulation of non-parametric estimation in which I need to minimize (maximize) a certain log-likelihood function. The function and commands I am using are:

LL=@(Bbin) (-1)*(sum(Yc'*log(proby*(fxnc./fxna))+(ones(M,1)-Yc)'*... log(ones(M,1)-proby*(fxnc./fxna))));

Bbin = fminsearch(LL,[-1;-1;-1])

As you can see, I am trying to find the Bbin that minimizes LL, but Bbin does not directly enter into my function "LL". Instead, it is used to calculate variables "fxnc" and "fxna". Those variables are calculated in such a way that I cannot add their calculation into the LL function (the code for how they are calculated is pasted below, but I don't believe it is relevant and may just add extra confusion). As you might guess, when trying to minimize LL, Matlab treats fxnc and fxna as an exogenous matrix, and therefore changing Bbin does not effect the LL function. Naturally, it spits out my original values for Bbin [-1;-1;-1]. How can I tell Matlab that varying Bbin changes fxnc and fxna?

Thanks!

Code for how fxnc and fxna are calculated:

%Calculate conditional means for Yc M=sum(Y); index=Xc*Bbin; si=std(index); fxnc=zeros(M,1); h=1.06*si*M^(-0.2); d=zeros(M,1); for j=1:M d(:,1)=(index(j,1)-index(:,1))./h; fx1c=normpdf(d); fxnc(j,1)=mean(fx1c/h); end

%Calculate conditional means for Y fxn=zeros(N,1); index=X*Bbin; si=std(index); h=1.06*si*N^(-0.2); d=zeros(N,1); for j=1:N d(:,1)=(index(j,1)-index(:,1))./h; fx1=normpdf(d); fxn(j,1)=mean(fx1/h); end

[Star Strider]

I can't figure out what you're doing, and I can't simulate your code to calculate means for Yc and Y, and LL. (Even when I plug in random numbers, it crashes. I can't get your d(:,1) assignment to work no matter what I try.) Since I can't test my ideas, this is a comment rather than an Answer.

It would help if you formatted your code, and provided some descriptions of the magnitudes and sizes of your variable vectors and matrices.

As a possible solution to your problem, I suggest you consider a for loop that defines Yc, Y, LL, and the rest, and iterates until you're happy with the results. You'll need to initialize outside the loop with guesses for the first pass.

[Jason]

Apologies, I should have included the entire code in my original post. Here it is below: ("if true" at the top and "end" at the bottom should be removed)

if true
clear all
N=200;
x1=random('unif',-3,-1,[N,1]);
x2=random('unif',-1,1,[N,1]);
x3=random('unif',2,4,[N,1]);
u=random('normal',0,1,[N,1]);
X=[x1 x2 x3];
B=[0.75;2;0.5];
Yn=X*B+2*u;
%Make Y Binary
Y=zeros(N,1);
for i=1:N;
  if Yn(i)>=0;
      Y(i)=1;
  else
      Y(i)=0;
  end
end
%Prob y=1
proby=sum(Y)/N;
%Create censored Y, where only Y=1 values exist
cen=[Y X];
cen(cen(:,1)==0,:)=[];
Yc=cen(:,1);
X1c=cen(:,2);
X2c=cen(:,3);
X3c=cen(:,4);
Xc=[X1c X2c X3c];
Bbin=[1;3;1];
%Calculate conditional means for Yc
M=sum(Y);
index=Xc*Bbin;
si=std(index);
fxnc=zeros(M,1);
h=1.06*si*M^(-0.2);
d=zeros(M,1);
for j=1:M
  d(:,1)=(index(j,1)-index(:,1))./h;
  fx1c=normpdf(d);
  fxnc(j,1)=mean(fx1c/h);
end
%Calculate conditional means for Y
fxn=zeros(N,1);
index=X*Bbin;
si=std(index);
h=1.06*si*N^(-0.2);
d=zeros(N,1);
for j=1:N
  d(:,1)=(index(j,1)-index(:,1))./h;
  fx1=normpdf(d);
  fxn(j,1)=mean(fx1/h);
end
%To evaluate, must eliminate points where Y=0
cen2=[Y fxn];
cen2(cen2(:,1)==0,:)=[];
fxna=cen2(:,2);
LL=@(Bbin) (-1)*(sum(Yc'*log(proby*(fxnc./fxna))+(ones(M,1)-Yc)'*...
  log(ones(M,1)-proby*(fxnc./fxna))));
Bbin = fminsearch(LL,[-1;-1;-1])

end

[Jan]

If "if true" should be removed, you could be so kind and remove by your own by editing the comment.

[Star Strider]

@Jason — I had time to experiment with your routine. Thank you for formatting it and providing the complete code.

One significant problem I discovered in looking through your code is that in your LL() function, Bbin doesn't actually appear anywhere as a parameter to be optimized! That's why it always returns the initial ‘guess’ values for it.

My next question: what precisely are you optimizing for? How do you define Bbin in terms of the rest of LL()? (NOTE that Bbin is a [3 x 1] vector and LL() produces a scalar.)

Also in LL(), note that the expression:

(ones(M,1)-Yc)'

is and will always be identically = 0 since Yc is a column vector of ones. So the second expression in LL()

(ones(M,1)-Yc)'*log(ones(M,1)-proby*(fxnc./fxna))

will therefore always be identically = 0 as well.

I'll keep this open for a bit while you sort out LL() and post a correction to it. I can't figure out what you're doing, so I can't correct it myself.

[Jason]

@Jan Simon- I'm new to this message board, but for some reason when I want to add code, "if true" appears. If I remove it, the code loses its formatting and is illegible. I'm not sure how to fix this or what the correct way to add code is, so I apologize for having to leave it there.

@Star Strider- You are correct when you point that ones(M,1)-Yc=0. I may have an error in how I'm defining my function, and will have to fix that- I'm not sure yet what the correct function should be. In the meantime, perhaps I can try to explain your first point. I understand that Bbin does not appear in the function, and therefore, of course, matlab cannot minimize LL by choosing Bbin. What I need matlab to do is minimize LL by indirectly choosing Bbin. Bbin effect fxna and fxnc, both of which appear in LL. I'm not sure how to code this though. If you (or anyone else) know of a way to do this, I can correct the LL function later on.

Thanks again for your help!

[Star Strider]

What do you want to do and how do you want Bbin to be optimised? What's your objective?

I'm lost.

[Jason]

Bbin should be chosen such that LL is minimized.

i.e. min (Bbin) LL(fxna(Bbin),fxnc(Bbin))

Answer 1

[Babak]

I wouldn't use an inline function like @(Bbin)... for this case.

I would do it this way: Write a function like this: (untested code - just a template)

 function LL_value = LL_function(Bbin)
 ...
 fxnc = ...
 fxna = ...
 LL_function = ...
 end

then in another script .m file, wirte routines that calls the LL_function recursively, like this:

 %initializing your parameters
 ...
 Bbin_initial = [-1;-1;-1]

 LL_value = LL_function(Bbin)

 Bbin_new = fminsearch(LL_function,Bbin_initial)

If this doesn't work fine, you can write your own minimization routine like this:

 %initializing your parameters
 ...
 Bbin_initial = [-1;-1;-1]

 while abs(Bbin_new - Bbin_old) > 0.001
   LL_value = LL_function(Bbin)

   Bbin_new = YourMinimizationRoutine(LL_function,Bbin_initial)
 end

Where you can use a newton raphson method to write the YourMinimizationRoutine function.

[Jason]

Thanks for your reply.I tried the first solution and it again just returned the initial values.

Regarding the second solution, I'm not sure how to apply my own minimization routine since there are three parameters to choose. I'm not familiar with the newton raphson routine (I'm a novice at Matlab). I've created minimization routines in the past for one parameter by creating a loop and adjusting the parameter value, but I'm not sure how to do this with multiple parameters in an efficient manner.

[Babak]

for three dimension, you can also write routines like a 1 dim minimizor. Try doing something like this:

 function Fvalue = func(vect)
    x = vect(1); y = vect(2); z = vect(3);
    Fvalue = x^2+y^2+z^2 % or the function you want to minimize
 end

then write your minimization script like MATLAB's minimizors, like,

 IC = [1,2,3];  %initial condition
 value_i = func(IC)
 delx_i=0.1;dely_i=0.1;delz_i=0.1;
 while abs(error)>1e-5 && i<10000

 % compute the value of the function as you deviate from x by delx_i
 x0 = IC+delx_i*[1,0,0];
 val_dx = func(x0)
 % compute the value of the function as you deviate from y by dely_i
 y0 = IC+dely_i*[0,1,0];
 val_dy = func(y0)
 % compute the value of the function as you deviate from z by delz_i
 z0 = IC+delz_i*[0,0,1];
 val_dz = func(z0)

 % Next, compare the new values val_dx,val_dy,val_dz,Value_i
and choose the direction at which the function value has decreased the most.  This direction is in the form of [a,b,c] like [3,-4,-12]/13. Then compute the value for the new input to the func you just found like IC+[a,b,c] and repeat until you either converge to a local/global minimum point or exceed a specific number of iterations,...

 end