[Solved] Fitting exponential decay function

329 views (last 30 days)
Grant Huckels
Grant Huckels on 25 Jul 2019
Edited: Image Analyst on 3 Feb 2020
I'm trying to fit an exponential decay to a dataset of x and y values (3001 each). Using other software I was able to calculate a k_off around 0.02 however using the fittype and fit to replicate this in MATLAB I get the following results:
Code:
s1 = sprintf('%f*exp(-koff*', y_equil); %(For y_equil = 0.148356)
s2 = 'x)+plateau'
eq_string = strcat(s1, s2);
f = fittype(eq_string);
f1 = fit(x,y, f)
plot(f1,x,y)
xlabel('Time (s)')
ylabel('Concentration')
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Results:
f1 =
General model:
f1(x) = 0.148356*exp(-koff*x)+plateau
Coefficients (with 95% confidence bounds):
koff = 5.888 (4.858, 6.917)
plateau = 0.1291 (0.1289, 0.1293)
This fit finds the correct value for the plateau but not for the k_off. For some reason the fit is horrible.
How do I correct this?
----------------------------------------------------------------
[Solved]: - I simply modified his code:
Credit to Jing Chen
https://www.mathworks.com/matlabcentral/fileexchange/37321-fitexponential-m
Corrected Code:
yBase = y_equil;
y = y - y_equil;
fh_objective = @(param) norm(param(2)+(param(3)-param(2))*exp(-param(1)*(x-x(1))) - y, 2);
initGuess(1) = -(y(2)-y(1))/(x(2)-x(1))/(y(1)-y(end));
initGuess(2) = y(end);
initGuess(3) = y(1);
param = fminsearch(fh_objective,initGuess);
k2 = param(1);
yInf = param(2) + yBase;
y0 = param(3) + yBase;
yFit2 = yInf + (y0-yInf) * exp(-k2*(x-x(1)));

Answers (2)

Image Analyst
Image Analyst on 26 Jul 2019
You can use fitnlm(). See attached demo. I can't use your actual data because you forgot to attach it.
00_Screenshot.png
  3 Comments
Walter Roberson
Walter Roberson on 3 Feb 2020
fit() with model 'exp' will automatically do the log transform internally.
However it you have a linear constant in the equation then you cannot use the log transform, except to get approximate starting values for a more accurate fit (assuming that none of the values are negative.)
Image Analyst
Image Analyst on 3 Feb 2020
Edited: Image Analyst on 3 Feb 2020
Knut, I'm not sure the coefficients would be the same if you did it that way as opposed to using the non-linear fit. Might - I'm just not sure.
It's easy to modify my code to get rid of the offset if you don't want it. Just change the formula and pass in only 2 values in the beta0 vector.
% Uses fitnlm() to fit a non-linear model (an exponential decay curve, Y = a * exp(-b*x)) 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 coefficients and 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.5 * 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.
% Note: it doesn't matter if X and Y are row vectors or column vectors since we use (:) to get them into column vectors for the table.
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));
% Guess values to start with. Just make your best guess.
aGuessed = 11 % Arbitrary sample values I picked.
bGuessed = 0.5
beta0 = [aGuessed, bGuessed]; % 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 = 30;
% Place formula text roughly in the middle of the plot.
formulaString = sprintf('Y = %.3f * exp(-%.3f * X)', coefficients(1), coefficients(2))
xl = xlim;
yl = ylim;
xt = xl(1) + abs(xl(2)-xl(1)) * 0.325;
yt = yl(1) + abs(yl(2)-yl(1)) * 0.59;
text(xt, yt, formulaString, 'FontSize', 25, 'FontWeight', 'bold');
% 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')
0001 Screenshot.png

Sign in to comment.


dpb
dpb on 26 Jul 2019
Edited: dpb on 26 Jul 2019
You forced the model to have f(0) = 0.148356 + 0.1291 = 0.2775.
It takes a really large decay constant to make up for that.
The initial value looks like it should be roughly that value of 0.148 so the multiplicating factor should be (0.148-plateau) if it is intended to fix that initial value rather than estimate three parameters instead of only two.
Try
y_equil=0.148356;
s=sprintf('(%f-plateau)*exp(-koff*x)+plateau', y_equil);
f = fittype(s)
f =
General model:
f(koff,plateau,x) = (0.148356-plateau)*exp(-koff*x)+plateau
>>
and see if don't have better luck
  1 Comment
Grant Huckels
Grant Huckels on 26 Jul 2019
Thank you for your time and effort! The k that I obtained using your code was around 0.8, and ~0.6 for the 'Robust' on.
However I found an alternative way thanks to Jing Chen's code. I left the updated material for future users that may encounter the same issues I did.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!