These things are always easier than you think. Well, except when they are just as hard as you think. ;-)
First, I'll create x as a COLUMN vector. That will avoid problems later, since you were forced to constantly transpose it. Just start out right from the beginning!
x = (0:100)';
Next, WHY IN THE NAME OF GOD AND LITTLE GREEN APPLES DID YOU USE RAND??????????????????
Do you understand that rand creates random numbers in the interval [0,1]? There will NEVER BE NEGATIVE NUMBERS COMING FROM RAND. USE RANDN INSTEAD!
y = exp(-0.1*x) + 0.1*rand(101,1); % THIS IS JUST BAD!!!!!
plot(x,y)
Next, plot it. Look at what you have. Note that here, the negative exponential is pretty much garbage noise for most of the curve, so the only useful signal comes from the initial part. But, ok. You still have enough data to get a fit out of it.
The model you have chosen is of the general form
y = a*exp(-k*x) + b
You MIGHT expect to find a==1, b == 0, and k = 0.1. Except that you used rand! And that means I will expect to see approximately b = 0.05.
Of course, those parameters will be estimated from somewhat poor data, but we will still get a valid result. Just for kicks, I'l throw the curve fitting toolbox at it first.
mdl = fittype('a*exp(-k*x) + b','indep','x')
fittedmdl = fit(x,y,mdl)
Even without starting values provided, fit was happy. We got something reasonable. Note that the estimate of b would be expected to be roughly 0.05. NOT zero. This is because the mean of a random number with a uniform distribution on the interval [0.0.1] is 0.05.
plot(fittedmdl,x,y)
I would not expect anything better. Now to look at your use of fminsearch for the fit. It appears you are trying to use a partitioned least squares, where you separate the conditionally linear parameters from the linear parameters. But you don't understand how to write that. (You could just download my fminspleas from the file exchange.)
k0 = 0.25;
[kest,abest]= pleasFit(k0, x, y)
As you can see, it returns the same estimates as did fit. Again, the estimate of b will be screwy because of your use of rand. The estimate I would espect to see for this model, with this data is ....
y = 1*exp(-0.1*x) + 0.05
function [kest,abest] = pleasFit(k0, x, y)
% Partitioned least squares estimation for the model
% y = a*exp(-k*x) + b
%
% k will be provided as an initial value, x and y are column vectors
linearest = @(x,y,k) [exp(-k*x),ones(size(y))]\y;
obj = @(k) norm([exp(-k*x),ones(size(y))]*linearest(x,y,k) - y);
options = optimset('PlotFcns',@optimplotfval,'Display','iter');
[kest,err,exitflag] = fminsearch(obj,k0,options);
abest = linearest(x,y,kest);
end
In the end, there were two MAJOR issues. First, WHY DID YOU USE RAND? RAND generates UNIFORM numbers in the interval [0,1]. I think too often, people just assume a random number is a random number, and do not recognize the difference between them. The result is you will get confusing results, because the additive random part does not have a mean of zero when you used rand.
Next, your code for the fit using fminsearch was completely wrong. It needs to estimate the parameters a and b INSIDE the objective. So for every value of k that fminsearch tries, it needs to estimate a NEW set of estimates for a and b. This was your major problem of course.
A partitioned nonlinear least squares splits the problem into two sets of unknowns. The nonlinear parameter (k), and then two conditionally linear parameters, here a and b. If you knew the value of k, you could estimate a and b. In fact, you could then just use polyfit to estimate a and b. (I probably should have done so.) But the linear part of the fit needs to happen INSIDE the objective function.
