user1018331
user1018331

Reputation: 143

vectorized function fminsearch

Thank you very much,first. The first code runs well:
function f=malibu(m1,m2,m3,k,l,t) f=(m1.*k+(m1+m2).*l).*exp(-m3.*t); end

function loglik= modelmalibu(p)
global k l t x n m2 m3;
f =malibu(p,m2,m3,k,l,t);
if f==0;
    loglik0=0;
else
loglik0=(x.*log(f)+(n-x).*log(1-f));%minus likelihood
end
loglik=sum(-loglik0);
end

clear all;
global n t x k l m2 m3;
m1=0.1;
m2=0.2;
m3=0.3;
t=[1 3 6 9 12 18]';
k=[1 1 2 3 3 4]';
l=[0 0 1 3 4 5]';
y=meltem(m1,m2,m3,k,l,t);
n=100;%trial
x=y.*n;%correct replies
pstart=0.3;
[p1,modelvalue]=fminsearch(@modelmalibu,pstart);

But the similar code for more variable gives error out.

function w=anemon(m1,m2,m3,X,Y,k,l)
w=(m1.*k+(m1+m2).*l)+X.*exp(-m3.*Y);
end

function loglik= modelanemon(p)
global n x m2 m3 X Y k l ;
f =anemon(p,m2,m3,X,Y,k,l);
if f==0;
    loglik0=0;
else
loglik0=(x*log(f)+(n-x)*log(1-f));%minus likelihood
end
loglik=sum(-loglik0);
end

clear;
global n x Ydata kdata ldata m1 m2 m3;
%parameters
m1=0.002;
m2=0.0001;
m3=7;
%given data
Xdata=[1 3 6 9 10 12]';
Ydata=[11 13 41 81 121 181]';
kdata=[1 1 2 4 5 4]';
ldata=[1 1 3 3 4 5]';
y=anemon(m1,m2,m3,Xdata,Ydata,kdata,ldata);
n=10;
x=y.*n;
pstart=2;
[pbest,modelvalue]=fminsearch(@modelanemon,pstart);

I've actually tried to use your advises, but if I would write an inequality instead of f==0, the first code fall down as well.

Upvotes: 0

Views: 990

Answers (1)

bdecaf
bdecaf

Reputation: 4732

I think you are mixing up two concepts.

vectorization: does refer in how to write your code so it can use some of the acceleration functions in your CPU. it has nothing to do with fminsearch. See: http://en.wikipedia.org/wiki/Vectorization_(parallel_computing)

I suppose you want to write your function such that it accepts a vector as input. Easiest way is to just use a function handle like this:

fh = @(x) my_complicated_function(const1, const2, x(1), x(2), x(3) )

In this case my_complicated_function has 5 inputs, and you take the first 2 constant and input a 3 dim vector for the other 3. Fminsearch will work with that.

You would call

 x_opt = fminsearch(fh, [1,2,3])

Besides some tips for the code:

  • Don't use == for comparison of numbers - go for eg. abs(x1-x2)<0.1 instead
  • interpanemon looks very strange - it doesn't do what one would call interpolation, and if you look - in every iteration p is recalculated so only the last calculation takes effect. As it looks it should output a constant value - no use optimizing that.
  • The use of a precalculated Z is probably not what you want. It already determines your optimum. If you use linear optimization the optimum must be already in Z - no need doing interpolation.
  • The invocation in your case might look like:

    min_p = fminsearch(@(x) interpanemon(5,Ydata,x,m2,m3,Z,X,Y,kdata,ldata) ,1)

Overall it looks very unusual - it could really help if you would explain the ideas WHY you choose to do it like this. Also it might help if you restate the problem in a simpler form.

Upvotes: 1

Related Questions