How to do exponential curve fitting like y=a*exp(b*x)+c

Hi guys,

I have a set of data x and y, which is given below: x=[10 12.5 15 17.5 20 22.5 25 27.5 30 32.5 35 37.5 40 42.5 45 47.5 50]; y=[62.1 77.3 92.5 104 112.9 121.9 125 129.4 134 138.2 142.3 143.2 144.6 147.2 147.8 149.1 150.9];

I'd like to to have a curve fitting like y=a*exp(b*x)+c. I tried to use cftool box (custom equation). However, it didn't work well. I am wandering if someone could help me with this.

Thanks

[Arturo Gonzalez]

Per this answer, you can do it as follows:

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  =  -1 + 5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);

% calculate integrals
iy1 = cumtrapz(x, y);
iy2 = cumtrapz(x, iy1);
iy3 = cumtrapz(x, iy2);

% get exponentials lambdas
Y = [iy1, iy2, iy3, x.^3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;

lambdas = eig([A(1), A(2), A(3); 1, 0, 0; 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -1.9997
%   0.5000

% get exponentials multipliers
X = [ones(size(x)), exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%  -0.9996
%   4.0043
%   1.9955
%   4.9999

Answer 1

[Star Strider]

Try this:

x=[10 12.5 15 17.5 20 22.5 25 27.5 30 32.5 35 37.5 40 42.5 45 47.5 50];
y=[62.1 77.3 92.5 104 112.9 121.9 125 129.4 134 138.2 142.3 143.2 144.6 147.2 147.8 149.1 150.9];

f = @(b,x) b(1).*exp(b(2).*x)+b(3);                                     % Objective Function
B = fminsearch(@(b) norm(y - f(b,x)), [-200; -1; 100])                  % Estimate Parameters

figure
plot(x, y, 'pg')
hold on
plot(x, f(B,x), '-r')
hold off
grid
xlabel('x')
ylabel('f(x)')
text(27, 105, sprintf('f(x) = %.1f\\cdote^{%.3f\\cdotx}%+.1f', B))

[MCC]

Hi Star Strider,

Thank you very much for your help. The code works well. Do you know how to determine the estimated parameters like [-200; -1; 100]?

Thanks

[Star Strider]

My pleasure.

Nonlinear parameter estimation is sensitive to the initial parameter estimates.

In this instance, I looked at the plotted data and guessed. It seems to have an upper asymptote of about 150, so ‘y=0’ would likely be about -150, and the exponent would be negative. My guess turned out to be accurate enough for fminsearch to converge on a reasonable set of estimated parameters.

In my posted code, I used initial parameter estimates that were close to the estimates fminsearch converged upon, using my previous guess as initial estimates.

[MCC]

Hi Star Strider,

Thank very much again. The curve fitting result is very sensitive to the initial value. I have been bothered by this problem for two days.

[Star Strider]

As always, my pleasure.

It is necessary to look at the data, and then estimate parameters from them. My initial guesses were:

[-100 -1 100]

Those worked. I posted initial estimates that are closer to the final estimates.

If my Answer helped you solve your problem, please Accept it!

[Gamra Samha]

Hi,

That is a very useful solution for plotting any exponential equaiton with a constant. I was wondering, would it then be possible to extract the goodness of fit such as rsquare, sse, and rmse values from this curve?

Thanks

Gamra

[Star Strider]

Some of them are available by using the fitnlm function. Others you would have to calculate yourself.

[Gamra Samha]

ok, will look into that. Thanks!

[Star Strider]

My pleasure.

[Alex Sha]

Guess of initial-start values is always a hard job for novices, in some cases, even for professional researchers, it is not easy work. For example, if the above fitting equation becomes form "y=b1*exp(b2*x)+b3" to "y=b1*exp(b2*x)+b3+b4*exp(b5/x)", it is almost impossible to get correct or near-correct initial-start values by manual, in this case, applying global optimization algorithms is maybe the only approach.

[Fatma Keskin]

I applied this to my own data but obtained a linear graph. I also changed the estimated parameters. I am wondering what the problem is actually.

[Kate Leary]

How would I calculate the confidence intervals for this fit? I was able to use the code to fit an exponential function to my data but I'd like to add confidence intervals.

[Star Strider]

@Kate Leary —

Use the fitnlm function —

x=[10 12.5 15 17.5 20 22.5 25 27.5 30 32.5 35 37.5 40 42.5 45 47.5 50];
y=[62.1 77.3 92.5 104 112.9 121.9 125 129.4 134 138.2 142.3 143.2 144.6 147.2 147.8 149.1 150.9];

f = @(b,x) b(1).*exp(b(2).*x)+b(3);                                     % Objective Function
% B = fminsearch(@(b) norm(y - f(b,x)), [-200; -1; 100])                  % Estimate Parameters
mdl = fitnlm(x(:), y(:), f, [-200; -0.1; 100])

mdl = 
Nonlinear regression model:
    y ~ b1*exp(b2*x) + b3

Estimated Coefficients:
          Estimate        SE         tStat       pValue  
          _________    _________    _______    __________

    b1      -203.88       4.5292    -45.013    1.5052e-16
    b2    -0.079833    0.0022031    -36.237    3.0597e-15
    b3       153.68      0.77353     198.67    1.4763e-25


Number of observations: 17, Error degrees of freedom: 14
Root Mean Squared Error: 1.01
R-Squared: 0.999,  Adjusted R-Squared 0.999
F-statistic vs. constant model: 5.57e+03, p-value = 4.9e-21

B = mdl.Coefficients.Estimate;
Bci = coefCI(mdl);
CoefficientMatrix = array2table([Bci(:,1) B, Bci(:,2)], 'VariableNames',{'-95% CI','Beta','+95% CI'})

CoefficientMatrix = 3×3 table
     -95% CI       Beta        +95% CI 
    _________    _________    _________

      -213.59      -203.88      -194.16
    -0.084558    -0.079833    -0.075108
       152.02       153.68       155.34

xv = linspace(min(x), max(x));
[yv,yci] = predict(mdl,xv(:));

figure
hp{1} = plot(x, y, 'pg', 'DisplayName', 'Data');
hold on
hp{2} = plot(xv, yv, '-r', 'DisplayName', 'Regression');
hp{3} = plot(xv, yci, ':r', 'DisplayName', '95% CI');
hold off
grid
xlabel('x')
ylabel('f(x)')
legend([hp{1};hp{2};hp{3}(1)], 'Location','best')
text(27, 105, sprintf('f(x) = %.1f\\cdote^{%.3f\\cdotx}%+.1f', B))

.

[Kate Leary]

@Star Strider

Thank you thank you! Is there a way to calculate prediction bounds as well? I know of predint but can't figure out how to use is without having a cfit or sfit object. Again, thank you for your help!

[Star Strider]

My pleasure!

They’re actually plotted (as ‘95% CI’, just so narrow that it’s difficult to see them here. the predict function call in my code calculates them:

[yv,yci] = predict(mdl,xv(:));

and they are then plotted in:

hp{3} = plot(xv, yci, ':r', 'DisplayName', '95% CI');

If you have any further questions, please feel free to follow up here.

.