Clear Filters
Clear Filters

Plateau followed by one phase decay

4 views (last 30 days)
Good morning, I am trying to figure out how to compute tau constants from my data
My data could be be fitted by such plateau followed by one phase decay function:
Here is an example:
x = 0:0.5:20; % time in seconds
Y0 = -0.6; % signal baseline value
Plateau = -1; % singnal plateu after trigger/stimulus, maximum change from baseline
tau = 0.6; % exponenential decay constant
K = 1/tau; % rate constant in units reciprocal of the x-axis units
X0 = 5; % trigger time
y = Y0*(x<=X0)+(Plateau+(Y0-Plateau)*exp(-K*(x-X0))).*(x>X0);
figure;plot(x,y,'k');
Given example raw data below, could you please support me with fitting of the function above and paramters extraction?
x = 0:0.5:20;
y = [-0.137055262721364 -0.118841612584876 -0.274602636741299 -0.117324828772196 ...
-0.173528150754918 -0.280491919000118 -0.244300356226590 -0.367583069701879 ...
-0.423274105143034 -0.529129050767333 -0.774173830727337 -0.676677606159725 ...
-0.730062482232667 -0.863905715495076 -0.831675679632950 -0.987303352625066 ...
-0.949979744575626 -0.865710605996821 -0.901728879393798 -0.877082148456042 ...
-0.944693953430828 -1.07404346760035 -0.915521627715257 -0.901789963321291 ...
-0.955365771797851 -0.941530617721837 -0.945983148775748 -1.01735658137382 ...
-0.965635004813717 -1.06321643780048 -0.956807780654745 -1.09208906741553 ...
-1.04341265165344 -1.08982901817714 -1.07984413818039 -0.934740294823467 ...
-0.960591807908718 -1.03623550995537 -0.909687220130007 -1.09290177705358 ...
-1.01208835337351];
Best regards
  2 Comments
Alex Sha
Alex Sha on 29 Apr 2024
refer to the results below:
Sum Squared Error (SSE): 0.160604734392522
Root of Mean Square Error (RMSE): 0.0625874479725772
Correlation Coef. (R): 0.979784143362497
R-Square: 0.959976967584583
Parameter Best Estimate
--------- -------------
x0 2.86487766740637
y0 -0.18364073498805
plateau -1.00952817360124
k 0.393751846945349
Francesco
Francesco on 29 Apr 2024
Thanks for your kind answer, would it be possible to share the code for fitting the function?
Best regards.

Sign in to comment.

Accepted Answer

Francesco
Francesco on 3 May 2024
Thanks a lot got for the quick reply, I am grateful to you and the community for the fantastic support.

More Answers (2)

Mathieu NOE
Mathieu NOE on 29 Apr 2024
hello
a code based on the poor's man optimizer (fminsearch)
no special toolbox required
% data
x = 0:0.5:20;
y = [-0.137055262721364 -0.118841612584876 -0.274602636741299 -0.117324828772196 ...
-0.173528150754918 -0.280491919000118 -0.244300356226590 -0.367583069701879 ...
-0.423274105143034 -0.529129050767333 -0.774173830727337 -0.676677606159725 ...
-0.730062482232667 -0.863905715495076 -0.831675679632950 -0.987303352625066 ...
-0.949979744575626 -0.865710605996821 -0.901728879393798 -0.877082148456042 ...
-0.944693953430828 -1.07404346760035 -0.915521627715257 -0.901789963321291 ...
-0.955365771797851 -0.941530617721837 -0.945983148775748 -1.01735658137382 ...
-0.965635004813717 -1.06321643780048 -0.956807780654745 -1.09208906741553 ...
-1.04341265165344 -1.08982901817714 -1.07984413818039 -0.934740294823467 ...
-0.960591807908718 -1.03623550995537 -0.909687220130007 -1.09290177705358 ...
-1.01208835337351];
% Fit the data with following fit model
% y = a if x<=x0
% y = b + (a-b)*exp(-c*(x-x0)) if x>x0
[a,b,c,x0] = plateauFit(x,y)
a = -0.1836
b = -1.0095
c = 0.3937
x0 = 2.8649
% Plot results
figure(1);
xa=linspace(min(x),x0,10); ya = a*ones(size(xa));
xb=linspace(x0,max(x),50); yb = b + (a-b)*exp(-c*(xb-x0));
plot(x,y,'k*',xa,ya,'-r',xb,yb,'-b')
grid on; xlabel('X'); ylabel('Y'); title('Special Fit')
xfit = [xa xb];
yfit = [ya yb];
[xfit,ia,ic] = unique(xfit);
yfit = yfit(ia);
yfit = interp1(xfit,yfit,x);
R2 = my_R2_coeff(y,yfit); % correlation coefficient
figure(2);
plot(x,y,'k*',x,yfit,'-r')
grid on; xlabel('X'); ylabel('Y');
title(['Special Fit / R² = ' num2str(R2) ], 'FontSize', 15)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [a,b,c,x0] = plateauFit(x,y)
%BILINEARFIT Fit plateau followed by exp decaying curve to a set of x,y data.
% The fitting function is
% y = a if x<=x0
% y = b + (a-b)*exp(-c*(x-x0)) if x>x0
% x, y should be vectors of equal length.
if length(x)~=length(y)
disp('Error in bilinearFit(x,y): x,y must have equal length.')
return
elseif min(size(x))>1 || min(size(y))>1
disp('Error in bilinearFit(x,y): x,y must be vectors.')
return
end
%% Determine initial guesses for a0, a1, b0, b1
% sort the data first
[x,I] = sort(x);
y = y(I);
%IC guess
a_init=y(1);
b_init=y(end);
x0_init = x(round(numel(x)/10));
% some special treatment for c_init
% assumed model : y = m*exp(c*x) once we have a bit transformed y data
yabs = abs(y);
yabs = yabs -yabs(1);
yabs = yabs(end) -yabs;
yabs(yabs<=0) = eps;
P = polyfit(x, log(yabs), 1);
c_init = -P(1);
f0=[a_init;b_init;c_init;x0_init]; % initial guess
% Define function myFun=difference between measured and predicted y
function v=myFun(f)
a=f(1); b=f(2); c=f(3); x0=f(4);
v=zeros(length(x),1);
for i=1:length(x)
if x(i)<=x0
v(i)= a -y(i) ;
else
v(i)= (b + (a-b)*exp(-c*(x(i)-x0))) -y(i);
end
end
end
f = fminsearch(@(x) norm(myFun(x)),f0);
a=f(1);
b=f(2);
c=f(3);
x0=abs(f(4));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R2 = my_R2_coeff(data,data_fit)
% R2 correlation coefficient computation
% The total sum of squares
sum_of_squares = sum((data-mean(data)).^2);
% The sum of squares of residuals, also called the residual sum of squares:
sum_of_squares_of_residuals = sum((data-data_fit).^2);
% definition of the coefficient of correlation (R squared) is
R2 = 1 - sum_of_squares_of_residuals/sum_of_squares;
end

Francesco
Francesco on 2 May 2024
Thanks a lot for your kind support with this fit issue.
I have tried the code above that you kindly provided on similar data, see attached data.mat
Unfortunately, the functions do not seem to work well with these type of data.
Your support and community support is highly appreciated and very helpful.
Best regards.
  1 Comment
Mathieu NOE
Mathieu NOE on 2 May 2024
seems we hve very tiny signals here , we can see the quantization effects
now to make the code work I needed to add some smoothing first
% data
load('data.mat');
% smooth a bit the data
ys = smoothdata(y,'gaussian',40);
% ys = smoothdata(y,'movmedian',40);
% Fit the data with following fit model
% y = a if x<=x0
% y = b + (a-b)*exp(-c*(x-x0)) if x>x0
[a,b,c,x0] = plateauFit(x,ys)
a = 1.0013
b = 0.9738
c = 11.0066
x0 = 1.0922
% Plot results
figure(1);
xa=linspace(min(x),x0,10); ya = a*ones(size(xa));
xb=linspace(x0,max(x),50); yb = b + (a-b)*exp(-c*(xb-x0));
plot(x,y,'k*',xa,ya,'-r',xb,yb,'-b')
grid on; xlabel('X'); ylabel('Y'); title('Special Fit')
xfit = [xa xb];
yfit = [ya yb];
[xfit,ia,ic] = unique(xfit);
yfit = yfit(ia);
yfit = interp1(xfit,yfit,x);
R2 = my_R2_coeff(y,yfit); % correlation coefficient
figure(2);
plot(x,y,'k*',x,ys,'b',x,yfit,'-r')
grid on; xlabel('X'); ylabel('Y');
title(['Special Fit / R² = ' num2str(R2) ], 'FontSize', 15)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [a,b,c,x0] = plateauFit(x,y)
%BILINEARFIT Fit plateau followed by exp decaying curve to a set of x,y data.
% The fitting function is
% y = a if x<=x0
% y = b + (a-b)*exp(-c*(x-x0)) if x>x0
% x, y should be vectors of equal length.
if length(x)~=length(y)
disp('Error in bilinearFit(x,y): x,y must have equal length.')
return
elseif min(size(x))>1 || min(size(y))>1
disp('Error in bilinearFit(x,y): x,y must be vectors.')
return
end
%% Determine initial guesses for a0, a1, b0, b1
% sort the data first
[x,I] = sort(x);
y = y(I);
%IC guess
a_init=y(1);
b_init=y(end);
% some special treatment for x0_init
thresh = 0.5*(y(1)+y(end));
[v,ind0] = min(abs(y-thresh));
x0_init = 0.5*x(ind0(1));
% x0_init = x(round(numel(x)/10));
% some special treatment for c_init
% assumed model : y = m*exp(c*x) once we have a bit transformed y data
yabs = abs(y);
yabs = yabs -yabs(1);
yabs = yabs(end) -yabs;
yabs(yabs<=0) = eps;
P = polyfit(x, log(yabs), 1);
c_init = -P(1);
f0=[a_init;b_init;c_init;x0_init]; % initial guess
% Define function myFun=difference between measured and predicted y
function v=myFun(f)
a=f(1); b=f(2); c=f(3); x0=f(4);
v=zeros(length(x),1);
for i=1:length(x)
if x(i)<=x0
v(i)= a -y(i) ;
else
v(i)= (b + (a-b)*exp(-c*(x(i)-x0))) -y(i);
end
end
end
f = fminsearch(@(x) norm(myFun(x)),f0);
a=f(1);
b=f(2);
c=f(3);
x0=abs(f(4));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R2 = my_R2_coeff(data,data_fit)
% R2 correlation coefficient computation
% The total sum of squares
sum_of_squares = sum((data-mean(data)).^2);
% The sum of squares of residuals, also called the residual sum of squares:
sum_of_squares_of_residuals = sum((data-data_fit).^2);
% definition of the coefficient of correlation (R squared) is
R2 = 1 - sum_of_squares_of_residuals/sum_of_squares;
end

Sign in to comment.

Categories

Find more on Data Type Identification in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!