Plateau followed by one phase decay
3 views (last 30 days)
Show older comments
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
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
Accepted Answer
More Answers (2)
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)
% 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
0 Comments
Francesco
on 2 May 2024
1 Comment
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)
% 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
See Also
Categories
Find more on Data Type Identification in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!