Exponential decay, rate constant

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!

[Mathieu NOE]

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) ?

[Anna Baj]

Hello,

yes, thank you. I think the problem were the upper and lower limits. I had to restrict them a bit.

I was wondering in terms of finding overall the best fit are there any guidelines that one should follow or just goodness of fit is enough?

Once again, thank you so much!

Answer 1

[Image Analyst]

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') 

[Image Analyst]

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

[Anna Baj]

Hi @Image Analyst, thank you so much.

Apologies for my mistakes, I am still learning.

What you posted here is amazing and it helped me understanding a lot more.

It's also good to have the chemical rate as I will need it also in the future!

I am extremely grateful!

Answer 2

[Mathieu NOE]

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');

[Mathieu NOE]

this code would be more adapted to a decaying exponential with four unknowns

% y = a.*(1-exp(b.*(x-c))) + d;

% generate some dummy noisy data
x = 0:100; 
N = length(x);
tau = 10; % time scale of relaxation 
y = 10-4*(1 - exp((-x-1)/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))) + 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_noisy);
sol = fminsearch(obj_fun, [y_noisy(end)-y_noisy(1),0,0,y_noisy(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_noisy,'r',x,y_fit,'-ok');
legend('signal','signal+noise','exp fit');

[Anna Baj]

Hello, thank you so much. I will try to adapt it to my data! I tried with many so far: a*exp(bx)+c

It could be that's I messed something up. I will let you know!

[Mathieu NOE]

ok

I can try on my side with your data if you wish

[Anna Baj]

If this wouldn't be a problem, it would be great.

Thank you!

[Mathieu NOE]

hello

here you are :

[numericData, ~, ~,] = xlsread('Bleaching.xlsx');
x  = numericData(:,1);
y2fit = numericData(:,2);
% exponential fit method 
% code is giving good results with template equation : % y = a.*(1-exp(b.*(x-c)));

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)-y2fit);
sol = fminsearch(obj_fun, [y2fit(end)-y2fit(1),0,0,y2fit(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,y2fit,'r',x,y_fit,'-.k');
legend('data','exp fit');

values for a,b,c,d are :

a_sol = -32.8773

b_sol = -0.0070

c_sol = 0.0047

d_sol = 152.5871

[Anna Baj]

Thank you so much! This looks really amazing. Thank you so much.

I realized that what I'd been doing was totally wrong.

[Mathieu NOE]

you're welcome !