Exponential decay, rate constant
71 views (last 30 days)
Show older comments
Hello,
I have a small problem. I'm trying to extract the rate of bleaching from my dataset.
It's a very simple dataset, yet I cannot make it work.
Best so far worked "nomal exp" but wasn't near perfect.
I would be extremely grateful for help.
Thank you!
2 Comments
Mathieu NOE
on 5 Feb 2021
hello again
I wonder why you use the same K parameter in your function for both amplitude and decay
shoudn't it be f(x) = A * exp (- K*x) ?
Accepted Answer
Image Analyst
on 6 Feb 2021
Anna: I think you're using the wrong equation. Your data clearly has an offset of about 120 so you need to add a constant to your equation. If you do that, you get a great fit. Try this with fitnlm():
% Uses fitnlm() to fit a non-linear model (an exponential decay curve, Y = a * exp(-b*x) + c) 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;
data = readmatrix('bleaching.xlsx');
whos data
% Create the X coordinates from 0 to 20 every 0.5 units.
X = data(:, 1);
Y = data(:, 2);
%--------------------------------------------------------------------------------------------------------------------------------------
% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b.', 'MarkerSize', 8);
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) + c
% 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)) + b(3);
% Guess values to start with. Just make your best guess.
aGuessed = 30 % Arbitrary sample values I picked.
bGuessed = 0.005
cGuessed = 120
beta0 = [aGuessed, bGuessed, cGuessed]; % 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) + coefficients(3);
% 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(-%.7f * X) + %.3f', coefficients(1), coefficients(2), coefficients(3))
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')
2 Comments
Image Analyst
on 6 Feb 2021
For what it's worth, I'm also attaching a demo for fitting a chemical reaction to the "Michaelis-Menten function" for chemical rate description. I know your data doesn't look like this, but it's a common situation and it may help others.
% Uses fitnlm() to fit data to the "Michaelis-Menten function" for chemical rate description.
% Initialization steps:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clearvars;
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 24;
%-----------------------------------------------------------------------------------------------------------------------------------
% Create sample data. You would not do this if you have your own data, obviously.
% Read in training data from a csv file.
baseFileName = 'fitnlm_rate_equation_data.csv';
if isfile(baseFileName)
% The file exists, so read it in and use the data from it.
data = csvread(baseFileName);
X = data(:, 1);
Y = data(:, 2);
else
% If file does not exist, create sample data from a function defined at the end of this script.
[X, Y] = CreateSampleData();
end
%-----------------------------------------------------------------------------------------------------------------------------------
% Plot the training data -- let's see what we're going to fit.
hFig2 = figure('Name', 'Demo by Image Analyst');
plot(X, Y, 'b.-');
grid on;
xlabel('X (Time)', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
title('Data to fit exponential to', 'FontSize', fontSize);
drawnow;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in. Use (:) to reshape data into column vectors.
tableXY = table(X(:), Y(:));
%-----------------------------------------------------------------------------------------------------------------------------------
% Fit a vertically shifted version of the the Fluorescence Recovery Curve (Michaelis-Menten function).
% One common non-linear equation is the Michaelis-Menten function,
% which is used to describe the reaction rate as a function of substrate concentration.
% The Michaelis-Menten function is characterized by a steep rise that gradually flattens into a plateau.
% Initially, small increases in substrate concentration cause a substantial increase in reaction rate,
% but eventually, additional increases in concentration cause progressively smaller increases in reaction rate.
% Eventually reaction rate plateaus and does not respond to increases in concentration.
% The Michaelis-Menten function has this form:
% rate = Vmax*C / (K + C)
% where r is the growth rate (the dependent variable),
% C is the concentration (the independent variable),
% Vm is the maximum rate in the system, and
% K is the concentration at which the reaction rate is Vm/2.
% These last two parameters (Vm and K) are the Michaelis-Menten parameters to be fit.
% Rewritten in X and Y, you get:
% Y = b1 * X ./ (b2 + X)
% Ref: http://strata.uga.edu/8370/lecturenotes/nonlinearRegression.html
% https://www.graphpad.com/guides/prism/7/curve-fitting/index.htm?reg_classic_hyperbola.htm
%-----------------------------------------------------------------------------------------------------------------------------------
% First we'll try that exact equation. After that, we'll try it again with an offset added to it to see if we can improve it.
modelFunction = @(b, tbl) b(1) * tbl(:, 1) ./ (b(2) + tbl(:, 1));
beta0 = [3, 1]; % Guess values to start with. Just make your best guess.
% Once you get close and see what the coefficients are, you can set beta0 to those values and
% Make another run to get better coefficients.
% Now the next line is where the actual model computation is done.
model1 = fitnlm(tableXY, modelFunction, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
coefficients1 = model1.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients1(1) * X ./ (coefficients1(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('Fitting Noisy Data to the Rate Equation', 'FontSize', fontSize);
drawnow;
% Compute residuals
residuals1 = sum(abs(Y - yFitted))
%-----------------------------------------------------------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% 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')
xticks(0:ceil(max(X))); % Put up tick mark every 1
xl = xlim; % Get the limits of the axes.
yl = ylim;
% Put the equation up on the figure.
xt = xl(1) + 0.11 * (xl(2) - xl(1));
yt = yl(1) + 0.50 * (yl(2) - yl(1));
message1 = sprintf('Model Equation 1 (has no offset): Y = %.3f * X / (%.3f + X)\nRate constant 1 = %.3f\nResiduals 1 Sum = %.1f', coefficients1(1), coefficients1(2), coefficients1(2), residuals1);
text(xt, yt, message1, 'Color', 'r', 'FontSize', fontSize);
%-----------------------------------------------------------------------------------------------------------------------------------
% Now let's try adding an offset to the equation to see if we can improve it (lower the residuals).
% It will work best for us if we add an offset to that, or
% Y = b1 * X ./ (b2 + X) + b3
modelFunction = @(b, tbl) b(1) * tbl(:, 1) ./ (b(2) + tbl(:, 1)) + b(3);
beta0 = [2.6, 1, 0.5]; % Guess values to start with. Just make your best guess.
% Once you get close and see what the coefficients are, you can set beta0 to those values and
% Make another run to get better coefficients.
% Now the next line is where the actual model computation is done.
model2 = fitnlm(tableXY, modelFunction, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
coefficients2 = model2.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients2(1) * X ./ (coefficients2(2) + X) + coefficients2(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
darkGreen = [0, 0.5, 0];
plot(X, yFitted, '-', 'Color', darkGreen, 'LineWidth', 2);
grid on;
drawnow;
% Compute residuals
residuals2 = sum(abs(Y - yFitted))
% Display the equation for model 2.
message1 = sprintf('Model Equation 2 (has offset): Y = %.3f + %.3f * X / (%.3f + X)\nRate constant 2 = %.3f\nResiduals 2 Sum = %.1f', coefficients2(3), coefficients2(1), coefficients2(2), coefficients2(2), residuals2);
yt = yl(1) + 0.250 * (yl(2) - yl(1));
text(xt, yt, message1, 'Color', darkGreen, 'FontSize', fontSize);
% Put up a line at the assymptote (what it would be at x = infinity where the ratio goes to 1).
assymptote1 = coefficients1(1)
yline(assymptote1, 'Color', 'm', 'LineWidth', 2);
assymptote2 = coefficients2(1) + coefficients2(3)
yline(assymptote2, 'Color', 'g', 'LineWidth', 2);
% Put up a legend.
legend('Noisy data', 'Model 1', 'Model 2', 'Assymptote for Model 1', 'Assymptote for Model 2', 'Location', 'east');
% Enlarge the tick labels.
ax = gca;
ax.FontSize = fontSize;
%====================================================================================================================
% Function to create sample data that follows the rate equation. This function is defined at the end of this script.
% Sample usage call:
% [x, yNoisy] = CreateSampleData();
% plot(x, yNoisy, 'b.');
% grid on;
function [x, yNoisy] = CreateSampleData()
% Create x values:
% Make 100 points in the 0-10 range (to give more precision in the steep part before the knee):
x1 = linspace(0, 10, 100);
% And make 250 points in the 10-100 range:
x2 = linspace(10, 100, 250);
% Now combine to give an overall x:
x = sort([x1, x2], 'ascend');
% Rescale x to be in the 0-15 seconds range.
x = rescale(x, 0, 15);
% Create a perfect y with no noise that follows the equation:
coefficients = [2.6, 1, 0.5]; % Guess values to start with. Just make your best guess.
% Create y with regular equation:
y = coefficients(1) * x ./ (coefficients(2) + x) + coefficients(3);
% Alternate way: Create y with anonymous function:
% modelFunction = @(b, x) b(1) * x ./ (b(2) + x) + b(3);
% y = modelFunction(coefficients, x) + 0.2 * rand(1, length(x));
% Add noise to the y values.
yNoisy = y + 0.2 * rand(1, length(x));
end
More Answers (1)
Mathieu NOE
on 5 Feb 2021
hello
this is a plain matlab example - i am not using the curve fitting tool here
% generate some dummy noisy data
x = 0:100;
N = length(x);
tau = 10; % time scale of relaxation
x_asym = 4; % relative increase of the observable x t -> infinity
y = x_asym*(1 - exp(-x/tau)); % uncorrupted signal
sig = 0.1; % noise strength
pert1 = sig*randn(1,N);
y_noisy = y+pert1; % corrupted signal
% exponential fit method
% code is giving good results with template equation : % y = a.*(1-exp(b.*(x-c)));
f = @(a,b,c,x) a.*(1-exp(b.*(x-c)));
obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y_noisy);
sol = fminsearch(obj_fun, [y_noisy(end),0,0]);
a_sol = sol(1);
b_sol = sol(2);
c_sol = sol(3);
y_fit = f(a_sol, b_sol,c_sol, x);
figure
plot(x,y,'-+b',x,y_noisy,'r',x,y_fit,'-ok');
legend('signal','signal+noise','exp fit');
7 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!