Fitting to exp plus a constant

Hello, I am trying to fit some data to an exponential + constant function and have this approach that I came across:

% code is giving good results with template equation : % y = a.*(1-exp(b.*(x-c))) + d;
f = @(a,b,c,d,x) a.*(1-exp(b.*(x-c))) + d;
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y);
sol = fminsearch(obj_fun, [y(end)-y(1),0,0,y(1)]);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
d_sol = sol(4);
y_fit = f(a_sol, b_sol,c_sol ,d_sol, x);
figure
plot(x,y,'-+b',x,y_fit,'-r'); grid on;
legend('signal','aexp(b(x-c))+d');
title(['a=',num2str(a_sol,'%0.2f'), '    b=',num2str(b_sol),'    c=',num2str(c_sol),'    d=',num2str(d_sol)],'FontSize',14,'Color','b')

Whilst the fit is good, the value of "a" seems wrong (at the expense of "d")

I've tried another approach that I came across from Steven Lord (as i don't really understand the code above)

D=[x y]
% It would be better to parametrize the function as 
xmin=min(x);
fitfun = fittype( @(a,b,c,x) a*exp(b*(x-xmin))+c );

[fitted_curve,gof] = fit(D(:,1),D(:,2),fitfun,'StartPoint',[max(y) -2 min(y)])
plot(fitted_curve,'r-')

But I get this error:

Error using fit>iFit (line 362)
Inf computed by model function, fitting cannot continue.
Try using or tightening upper and lower bounds on coefficients.
Error in fit (line 117)
[fitobj, goodness, output, convmsg] = iFit( xdatain, ydatain, fittypeobj, ...
                                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in PlotAnalysis/FitExp1ContstantButtonPushed (line 1361)
[fitted_curve,gof] = fit(D(:,1),D(:,2),fitfun,'StartPoint',[max(y) -2 min(y)])
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Answer 1

[Matt J]

It looks like the fit was successful, but the model function is overparametrized, so there is no specific value for a (or c) that you can count on,

where . Whatever the decay rate b needed for the fit, there are infinite choices for a and c that will produce the same A (and thus the same curve).

[Jason]

thanks Matt. So looking at the curve, the constant value d is definetly below the last data point. This would then mean "a" would be positive . How can I force the "d" term to be less than the last data point?

[Jason]

Hmmm, I can't see my obvious mistake!

f = @(a,b,d,x) a.*(exp(-b.*x)) + d;
obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);
sol = fminsearch(obj_fun, [y(end)-y(1),0,0,y(1)]);
a_sol = sol(1);
b_sol = sol(2);
d_sol = sol(3);
y_fit = f(a_sol, b_sol ,d_sol, x);

[Matt J]

This is the model you show in your original code,

f = @(a,b,c,d,x) a.*(1-exp(b.*(x-c))) + d;

This model spans the same space of curves as a*exp(-b*c)+d, but with this parametrization, it is appropriate that "a" would be negative.

[Sam Chak]

Hi @Jason

To understand what occurred, you need to plot both models for the negative and positive signs of the function's amplitude (a).

One model exhibits exponential growth, while the other demonstrates exponential decay. The constant parameter "d" simply shifts the exponential curve upward or downward.

% domain
x   = 0:.01:100;

% parameters
a   = -2.58;
b   = -0.020469;
c   =  0.020776;
d   =  750.5648;    % bias

% models
z   = exp(b*(x - c));       % pure exponential model
y1  =  a*(1 - z) + 0;       % exponentially decaying
y2  = -a*(1 - z) + 0;       % exponentially growing

% plot result
figure
plot(x, [y1; y2]), grid on
xlabel('x')
ylabel('Amplitude')
legend('exponentially decaying', 'exponentially growing', 'location', 'east')

●

figure
plot(x, z), grid on
xlabel('x')
ylabel('Amplitude')
title('decaying exponential model: exp(- 0.020469*(x - c))')

[Matt J]

@Jason Just to put a finer point on it, there is no reason to think the value of "d" you're seeing is inappropriate in any way, or linked to the negative sign you're seeing on the "a" parameter. We know the sign of the "a" parameter is correct because otherwise, your convex data would end up being fitted by a concave curve. That's true independently of the specific form of the model equation that you're using.

[Jason]

Thankyou all.

I have decided this way works and it easier to read (I couldn't get the other apporach to work)

% 2nd approach
x1=x; y1=y;
fitfun = fittype( @(a,b,c,x) a*exp(-b*x)+c);
[fitted_curve,gof] = fit(x1,y1,fitfun,'Lower',[0,0,0])      % 'StartPoint',[max(y) -2 min(y)
hold(ax,'on');
p=plot(ax,fitted_curve,'r--',x1,y1,'y--');   % 
% p.DisplayName='aexp(-bx)+c'

But some odd issues came up. My raw data x & y is already ploted on a UIAxes - ax.

It seems you cannot do this

hold(ax,'on');
p=plot(ax,fitted_curve,'r--');   % 

Instead you have to include the raw data too

hold(ax,'on');
p=plot(ax,fitted_curve,'r--',x1,y1,'y--');   % 

So now I have 2 plots of the same data.

It also doesn't allow this to work

 p.DisplayName='aexp(-bx)+c'   % Add to legend.
 

Am I missing something to be able to just add the fitted curve alone to my current UIAxes?

[Jason]

Found the solution! (didn't know about the formula function!

p=plot(ax,x1,fitted_curve(x1),'y--'); %

f=formula(fitted_curve)

p.DisplayName=f

[Jason]

This is really interesting. i was about to write up my compeete soluton and have noticed that using this approach, sometimes the fits are way off, it seems a bit random.

Here is the code im using.

ax=app.UIAxes; clc;
[x1,y1,cl] = getDataFromGraph(app,ax,1);    % My function toe xtract x, y data from plot


% 2nd approach
fitfun = fittype( @(a,b,c,x) a*exp(-b*x)+c);
[fitted_curve,gof] = fit(x1,y1,fitfun,'Lower',[0,0,0])      % 'StartPoint',[max(y) -2 min(y)

hold(ax,'on');
p=plot(ax,x1,fitted_curve(x1),'y--');   % plot fit
f=formula(fitted_curve);
p.DisplayName=f;

% C=coeffvalues(fitted_curve)
% a=C(1); b=C(2); c=C(3);
% title(ax,['a=',num2str(a,'%0.2f'), '    b=',num2str(b),'    c=',num2str(c)],'FontSize',18,'Color','c')

N=coeffnames(fitted_curve)
V=coeffvalues(fitted_curve)

n=numel(N)
s=[];
for i=1:n
    snew=[N{i},'=',num2str(V(i)),'   '];
    s=[s snew]

end
title(ax,s,'FontSize',18,'Color','c')

Here are some a,b,c values from fits

a=2.5846 b=0.020469 c=747.9813 --> Good fit

a=108.4951 b=0.00020374 c=641.6853 --> Bad fit

(nothing changes so it seems hit and miss)

[Matt J]

You should not use a custom fittype model for the fit. This requires you to provide a very accurate initial guess of the parameters. Instead, use the built-in 'exp2' fittype, with bounds to constrain the second decay rate to zero.

Alternatively, fminspleas from the File Exchange would also be very robust,

https://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

a = 2.5846; b = 0.020469; c = 747.9813;  %hypothetical data
x=1:100;
y=a*exp(-b*x)+c; y=y+randn(size(y))*0.1;

[b,ac]=fminspleas( {@(b,x) exp(-b*x),1} ,0.005, x,y,0); %Do fit
a=ac(1); c=ac(2);

plot(x,y,'--bx',x, a*exp(-b*x)+c,'-r'); %Plot
title(sprintf('a=%.02f, b=%.04f,c=%.02f,',a,b,c))
legend Signal 'a*exp(-b*x)+c'

●

[Torsten]

Supply initial values for the fitting parameters. Then you should get back the same results.

Maybe "fit" uses randomized initial guesses.

[Jason]

Hi Matt, could I ask what the 0.005 is in this please, and where the constants a & c are

[b,ac]=fminspleas( {@(b,x) exp(-b*x),1} ,0.005, x,y,0); %Do fit

[Matt J]

I think you would get a better handle on that if you familiarized yourself with the fminspleas documentation. Briefly, though,

(1) 0.005 is the initial guess for b

(2) a and c are the coefficients applied to the exp(-b*x) and the 1 respectively in the input {@(b,x) exp(-b*x),1}. Their fitted values are returned in "ac". They do not require initial guess values.