MATLAB Answers

How do I find the coefficients of a function for a curve fit?

170 views (last 30 days)
I have the data:
X = [0; 0.5; ...]
Y = [0; 0.8; ...]
And the function to which I want to model the data:
F = a + exp ^ (- x * b)
How do I find "a" and "b" for a good curve fit?
I should not use cftool.

  7 Comments

Show 4 older comments
Martín Menéndez
Martín Menéndez on 22 Apr 2017
Thanks for you help an time
I removed terms to the function so it was not so long;
clc;
close all;
clear;
workspace;
format long g;
format compact;
fontSize = 20;
load ('X.txt')
load('Y.txt')
X=X';
Y=Y';
a=10;
b=6;
c=4;
d=2;
e=3;
Y= a-b*exp(-c*X)-d*exp(-e*X);
Y = Y + 0.05 * randn(1,length(Y));
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
tbl = table(X', Y');
modelfun = @(b,x) b(1)-b(2)*exp(-b(3)*x(:, 1))-b(4)*exp(-b(5)* x(:, 1));
beta0 = [10, 6, 4, 2, 3];
mdl = fitnlm(tbl, modelfun, beta0);
coefficients = mdl.Coefficients{:, 'Estimate'};
yFitted = coefficients(1) - coefficients(2)* exp(-coefficients(3)*X) - coefficients(4)* exp(-coefficients(5)*X);
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
Walter Roberson
Walter Roberson on 23 Apr 2017
> In nlinfit>LMfit (line 579)
In nlinfit (line 276)
In NonLinearModel/fitter (line 1123)
In classreg.regr.FitObject/doFit (line 94)
In NonLinearModel.fit (line 1430)
In fitnlm (line 94)
In guava (line 49)
Warning: Rank deficient, rank = 3, tol = 1.458001e-13.
> In nlinfit>LMfit (line 579)
In nlinfit (line 276)
In NonLinearModel/fitter (line 1123)
In classreg.regr.FitObject/doFit (line 94)
In NonLinearModel.fit (line 1430)
In fitnlm (line 94)
In guava (line 49)
Warning: Rank deficient, rank = 3, tol = 1.451491e-13.
> In nlinfit>LMfit (line 579)
In nlinfit (line 276)
In NonLinearModel/fitter (line 1123)
In classreg.regr.FitObject/doFit (line 94)
In NonLinearModel.fit (line 1430)
In fitnlm (line 94)
In guava (line 49)
Warning: Rank deficient, rank = 3, tol = 1.450838e-13.
Warning: Some columns of the Jacobian are effectively zero at the solution, indicating that the model is insensitive to some of its parameters. That may be
because those parameters are not present in the model, or otherwise do not affect the predicted values. It may also be due to numerical underflow in the model
function, which can sometimes be avoided by choosing better initial parameter values, or by rescaling or recentering. Parameter estimates may be unreliable.
> In nlinfit (line 373)
In NonLinearModel/fitter (line 1123)
In classreg.regr.FitObject/doFit (line 94)
In NonLinearModel.fit (line 1430)
In fitnlm (line 94)
In guava (line 49)

Sign in to comment.

Accepted Answer

Image Analyst
Image Analyst on 22 Apr 2017
If you don't have the Curve Fitting Toolbox but have the Statistics and Machine Learning Toolbox, you can use fitnlm():
% Uses fitnlm() to fit a non-linear model (an exponential decay curve) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% Create the X coordinates from 0 to 20 every 0.5 units.
X = 0 : 0.5 : 20;
% Define function that the X values obey.
a = 10 % Arbitrary sample values I picked.
b = 0.4
Y = a + exp(-X * b); % Get a vector. No noise in this Y yet.
% Add noise to Y.
Y = Y + 0.05 * randn(1, length(Y));
% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X', Y');
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) + exp(-b(2)*x(:, 1));
beta0 = [10, .4]; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) + exp(-coefficients(2)*X);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 25;
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

  0 Comments

Sign in to comment.

More Answers (2)

Walter Roberson
Walter Roberson on 22 Apr 2017
The best I find for the 7-coefficient model, using custom code that I have not released, is
[1.95345951298651754, -0.271402425481958642, 0.114129599568334961, 0.6048791673085121, 0.000138307518530737773, 1.61189969373566888, -0.0124593500324885736]
You can get much the same result using cftool if you adjust the bounds. For example for each positive number above, set the lower bound to 0 and upper to 2; and for each negative number above set the lower bound to -2 and the upper to 0.
I evaluated over a large range of coefficients; there are some other locations that are not bad, but this was the best I found. For example, at
[1.95618718619627652, -0.206949316792977123, 48925.8303936270968, 1.56148395793133599, 0.011815555675624435, 0.601652543952944985, -0.000134723232812764676]
the residue is only about 15% larger.

  1 Comment

Walter Roberson
Walter Roberson on 23 Apr 2017
For your reduced system with 5 coefficients,
Y= a-b*exp(-c*X)-d*exp(-e*X)
there are two minima that are equivalent, one near
[1.96086280008987 0.596272887303567 0.000128842293439593 1.50536232706941 0.0109031671832912]
and the other near
[1.96086281505995 1.50536227966151 0.0109031664979331 0.596272910300341 0.000128842291110631]
You can get close to the first of those in a fraction of a second using cftool custom equation and using the Fit Options to constrain each of the parameters between 0 and 2.

Sign in to comment.


Alex Sha
Alex Sha on 24 Feb 2020
The below should be the best solution:
Root of Mean Square Error (RMSE): 0.0255893453091926
Sum of Squared Residual: 0.047146650721423
Correlation Coef. (R): 0.999022475612635
R-Square: 0.998045906779199
Adjusted R-Square: 0.997989266395987
Determination Coef. (DC): 0.998045906779199
Chi-Square: 0.0562099503181627
F-Statistic: 5533.08505193596
Parameter Best Estimate
---------- -------------
a 1.95345951949563
b -0.271402458132563
c 0.114129545758927
d 1.61189973355979
e 0.0124593503150589
f 0.604879174590309
g -0.000138307512053813

  0 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!