Is this the correct optimization tool (lsqcurvefit,MultiStart)

I have checked the Optimization tool box to see which one would be good for my optimization problem, so I thought lsqcurvefit would be good for my problem. Upon making a model with data with known parameters and trying to fit it with a model. The parameters estimated are not the same as the known parameters (or converge to the known parameters) that I first created the noisy data. What do you think is wrong with this? I'm also using MultiStart to find different initial points also.

%Creating data with noise:
tu is a (225x3) data 
tdata is 225x1 =tu(:,1);
ptest=[0.2,0.1,0.25,0.05];
Tissue= odefnc_noise(ptest,tdata,tu);

%Fitting the data with noise
fun = @(p,tdatas) odefnc_test(p,tdatas,tu);
param0=[0.1 0.1 0.1 0.1];
problem = createOptimProblem('lsqcurvefit','objective',fun,...
'xdata',tdata,'ydata',Tissue,'x0',param0,'lb',[0 0 0 0],'ub',[1 1 1 1],...
'options',options);
[b,fval,exitflag,output,solutions]=run(ms,problem,2);

%Fitting Function
function Tissue= odefnc_test(p,texp,tu)
F=p(1);
fp=p(2);
fis=p(3);
PS=p(4);
init = [0 0];
[~,y]=ode45(@(t,x) myeqn(t,x,F,fp,fis,PS,tu), texp, init);
P=y(:,1);
I=y(:,2);
Tissue = (I+P);

function dx=myeqn(t,x,F,fp,fis,PS,tu)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1)= (F/fp)*(output-x(1))- (PS/fp)*(x(1)-x(2));
dx(2)=  (PS/fis) *(x(1)-x(2));
dx=[dx(1);dx(2)];
end
end

This will converge to different value even if param0=ptest. So is it the ode function that create a problem? If 'fun' is not an ode, so it is like fun=@(b) b(1)*xdata +b(2)*xdata^2, the solution would converge according to many other problems I found online.

[Matt J]

If the data is noisy, one can't expect the true parameters to be exactly recovered. There will be some error in the estimate due to the noise.

[Arbol]

So it will never or rarely converge even if the noise is gaussian noise?

[Walter Roberson]

Right, even with gaussian noise the recovered parameters will seldom exactly equal the original parameters.

[Walter Roberson]

Please post the code for odefnc_noise; also if you could post tu and tdata values for us to test with.

[Arbol]

This is my odefnc_noise:

function Tissue= odefnc_noise(p,texp,tu)
F=p(1);fp=p(2);fis=p(3);PS=p(4);
init = [0 0];

[~,y]=ode45(@(t,x) myeqn(t,x,F,fp,fis,PS,tu), texp, init);
P=y(:,1);
I=y(:,2);
p=max(P)-min(P);
i=max(I)-min(I);
b=mean(p+i);
noise= 0.01*b*randn(size(texp));
Tissue = (I+P)+noise;

function dx=myeqn(t,x,F,fp,fis,PS,tu)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1)= (F/fp)*(output-x(1))- (PS/fp)*(x(1)-x(2));
dx(2)=  (PS/fis) *(x(1)-x(2));
dx=[dx(1);dx(2)];
end
end

And I have attached my test.txt, tu= load test. And tdata = tu(:,1); And if you want to test the real data, instead of the simulated data from odefnc_noise, the data is from tu(:,3).

[Arbol]

@Walter Roberson, did you find any mistakes on my code? thanks

[Walter Roberson]

We do not have your options or your ms

You have

Tissue = (I'+P')+noise;

Here, I, P and noise are all column vectors. The (I'+P') would be a row vector. You would then be adding a row vector and a column vector. In R2016a and earlier that would be an error; in R2016b and later it is the equivalent of

bsxfun(@plus, (I'+P'), noise)

which would produce a 225 x 225 output. The 225 x 225 output would be a problem in later steps.

[Arbol]

Here are the options:

options = optimset('Display','iter','FinDiffRelStep',0.0015,'TolFun', 2e-15,...
'TolX',2e-15,'DiffMaxChange',1,...
    'FunValCheck','on',...
    'MaxFunEvals',300,...
    'MaxIter',300);

I have played around with 'FinDiffRelStp', and found that it is faster and converge (not flag=1) to a value rather than the default setting. I also relaxed my 'TolFun' and 'TolX'. Whenever I run this, I always get exitflag=0. Also, I have fixed my noise function to column+column in the previous code.

[Walter Roberson]

Please post your construction of ms which is probably a MultiStart object.

Also, to check, did you change

Tissue = (I'+P')+noise;

to

Tissue = (I+P)+noise;

?

[Arbol]

I have ms as:

ms=MultiStart;
ms.UseParallel='always';

I did check with different Tissue, they have the same result. The fit looks good but it doesnt have exitflag=1.

[Arbol]

Also, I did change Tissue = (I+P)+noise;

function Tissue= odefnc_noise(p,texp,tu)
F=p(1);fp=p(2);fis=p(3);PS=p(4);
init = [0 0];
[~,y]=ode45(@(t,x) myeqn(t,x,F,fp,fis,PS,tu), texp, init);
P=y(:,1);
I=y(:,2);
p=max(P)-min(P);
i=max(I)-min(I);
b=mean(p+i);
noise= 0.01*b*randn(size(texp));
Tissue = (I+P)+noise;
function dx=myeqn(t,x,F,fp,fis,PS,tu)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1)= (F/fp)*(output-x(1))- (PS/fp)*(x(1)-x(2));
dx(2)=  (PS/fis) *(x(1)-x(2));
dx=[dx(1);dx(2)];
end
end

[Walter Roberson]

Your tolerance is set to 2E-15, but your function values are coming out on the order of 50000 or 100000 in the areas being explored, so your tolerance is much lower than eps() of the function values. You will therefore be unable to converge and you will just exhaust your maximum function calls or maximum interations. Then the multistart process will notice that the calls all failed to converge and will refuse to run any more locations.

[Arbol]

So does this means that I have to increase my tolerance? If I do, I will flag 2 or 3 due to step size too small or have reached the tolerance level.

[Arbol]

I'm sorry, here try this new codes that I did:

%----------------Creating data with noise:------------------%
tu is a (225x3) data 
tdata is 225x1 =tu(:,1);
ptest=[0.2,0.1,0.25,0.05];
Tissue= odefnc_noise(ptest,tdata,tu);

%-------------Fitting the data with noise--------------%
options = optimset('Display','iter','FinDiffRelStep',0.0015,'TolFun', 2e-15,...
'TolX',2e-15,'DiffMaxChange',1,...
  'FunValCheck','on',...
  'MaxFunEvals',300,...
  'MaxIter',300);
ms=MultiStart;
ms.UseParallel='always';
fun = @(p,tdatas) odefnc_test(p,tdatas,tu);
param0=[0.1 0.1 0.1 0.1];
problem = createOptimProblem('lsqcurvefit','objective',fun,...
'xdata',tdata,'ydata',Tissue,'x0',param0,'lb',[0 0 0 0],'ub',[1 1 1 1],...
'options',options);
[b,fval,exitflag,output,solutions]=run(ms,problem,2);

%------Fitting Function--------%
function Tissue= odefnc_test(p,texp,tu)
F=p(1);
fp=p(2);
fis=p(3);
PS=p(4);
init = [0 0];
[~,y]=ode45(@(t,x) myeqn(t,x,F,fp,fis,PS,tu), texp, init);
P=y(:,1).*fp;
I=y(:,2).*fis;
Tissue = (I+P);
function dx=myeqn(t,x,F,fp,fis,PS,tu)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1)= (F/fp)*(output-x(1))- (PS/fp)*(x(1)-x(2));
dx(2)=  (PS/fis) *(x(1)-x(2));
dx=[dx(1);dx(2)];
end
end

%---------Data with Noise------%:

function Tissue= odefnc_noise(p,texp,tu)
F=p(1);fp=p(2);fis=p(3);PS=p(4);
init = [0 0];
[~,y]=ode45(@(t,x) myeqn(t,x,F,fp,fis,PS,tu), texp, init);
P=y(:,1).*fp;
I=y(:,2).*fis;
p=max(P)-min(P);
i=max(I)-min(I);
b=mean(p+i);
noise= 0.01*b*randn(size(texp));
Tissue = (I+P)+noise;
function dx=myeqn(t,x,F,fp,fis,PS,tu)
output=interp1(tu(:,1),tu(:,2),t); 
dx(1)= (F/fp)*(output-x(1))- (PS/fp)*(x(1)-x(2));
dx(2)=  (PS/fis) *(x(1)-x(2));
dx=[dx(1);dx(2)];
end
end

Just a minor change to P and I from the older code; changed to P=y(:,1)*fp and I=y(:,2)*fis respectively. I included the ms and options. Thank you for the follow up, I appreciate it.

[Walter Roberson]

I have not run with your modifications yet. With the version before them, when I changed the tolerance and step size tolerance to some that looked reasonable, nearly all of the runs end with "Change in x was less than the specified tolerance." which is exit 2; the remaining run was at a higher function value and ran out of iterations (exit 0).

I pushed up the number of starts to 20, but by random chance the final output was not better than one of the runs with only 2 multistarts -- so the starting position is important.

[Arbol]

Correct, the number of points doesn't make it any better. But the new correction code that I posted will make the fit better, but still doesn't give a good exitflag. Even when i raised the number of starts.