# Plateau followed by one phase decay

4 views (last 30 days)
Francesco on 29 Apr 2024
Answered: Francesco on 3 May 2024
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
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 on 29 Apr 2024
Thanks for your kind answer, would it be possible to share the code for fitting the function?
Best regards.

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.

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 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.
Best regards.
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
% 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

### Categories

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

R2019a

### Community Treasure Hunt

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

Start Hunting!