Reputation: 143
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
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:
==
for comparison of numbers - go for eg. abs(x1-x2)<0.1
insteadinterpanemon
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 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