Plateau followed by one phase decay

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]

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]

Thanks for your kind answer, would it be possible to share the code for fitting the function?

Best regards.

Answer 1

[Francesco]

Thanks a lot got for the quick reply, I am grateful to you and the community for the fantastic support.

Answer 2

[Mathieu NOE]

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

Answer 3

[Francesco]

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.

[Mathieu NOE]

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