I want to fit the curve with two linear curves and find the point where the fitted curve deviates from the linear relationship, sort of the tangent point.

I couldn't find an appropriate answer to my question. I have attached the matlab script which plots the results and also the clean data. i want to fit two linear curves and identify the slopes. I tried depicting in the screenshot. Additionally i want to know the point where the fit deviates from the curve.

Maybe easier way is to split the data where we see two different scales. Fit fast rise region and then slow rise region.

Please suggest! Appreciate your help.

[Bruno Luong]

Frankly it needs imagination to see those curvy data fitted to 2 straight lines.

[William Rose]

@Bruno Luong is right. I would not use two straight lines for these data.

Answer 1

[Image Analyst]

Find the point where the slope difference between the left and right side of the dividing point is greatest.

See my attached demo. Adapt it to use your data instead of my demo data.

[T]

Nice! Thanks @Image Analyst. I tried this. This works well for the fitting at linear x and y scale as you show above but I'm struggling to fit the data on a log-log scale with the attached data.

[William Rose]

This is a nice script from @Image Analyst!

For log-log data, just do the fit on logx=log(x) versus logy=log(y).

Note that the approach of @Image Analyst, like bilinearFit, does a least-squares fit of each segment. For bilinearFit, you said that you wanted to find the tangent line rather than a line through the data. Of course the lines found by @Image Analyst and by bilinearFit have a slope that is tangent to the segment somewhere, if you shift them a bit.

Why do you want to maximize the difference of the two fitted slopes? Is there a modeling reason, or an estimation argument, that says maximizing this difference leads to desired result?

If I were reviewing a manuscript or dissertation that used this approach, I would ask why the curve should be fit with two lines, instead of one, or three, or more. The fit gets better when you use more than one line segment, of course. Does the improvement in fit justify the extra parameters? The Akaike information criterion is a classic way to answer this question. See here.

[T]

Thanks @William Rose, you have been very helpful! I still dont have the result I would like, something like I show in the attached image.

This function identifies the point at ~(5,-4) but I would expect a point at ~(2.5,-4).

Really, I am not looking for a very strong reasoning for drawing two tangent curves at the change of slope on log-log plot. This is for a hypothesis I'm working on where eventually I would eventually ignore the second part of the curve (which is almost a flat line). So essentially I want to find the intersection point of two tangents of these two curves and while finding the slope of the first curve on log-log, ignore the second half of the curve. Sorry it might be difficult to understand but I am just trying something very basic! There wont be a review for this work :p, it's just a personal study :D Many thanks!

[Image Analyst]

@T You need to widen your search range. Look at my code:

% Assume the crossing point will be somewhere in the middle half of the points.
% If you go much more outward than that there are too few points to get a good line and the slopes of the lines will vary tremendously.
% Fit a line through the right and left parts and get the slopes.
% Keep the point where the slope difference is greatest.
numPoints = length(x);
index1 = round(0.25 * numPoints); % 25% of the way through.
index2 = round(0.75 * numPoints); % 75% of the way through.
% In other words, assume that we need at least 25 percent of the points to make a good estimate of the line.
% Obviously if we took only 2 or 3 points, then the slope could vary quite dramatically,

Try going from 10% to 90% instead of from 25% to 75%.

[William Rose]

@T,

[edit: change "there" to "therefore" near the end of my comments]

Sorry I went a bit academic in my previous comment. That is a mistake I sometimes make.

You write: "eventually I would eventually ignore the second part of the curve (which is almost a flat line). So essentially I want to find the intersection point of two tangents of these two curves and while finding the slope of the first curve on log-log, ignore the second half of the curve.". This is interesting and I understand it. The tricky part about doing this is deciding what part of the curve to ignore. You said "second half", but I assume you did not necessarily mean exactly half.

I suggest that you fit the data with two straight lines, and require the right hand line to be flat. Then see where the left line intersects the flat right line.

There's still a question, if you do the approach I suggest: Do you want the flat right line's y-axis intercept to be (A) at the y-asymptote of the curve, or do you want the flat right line's y-intercept to be (B) fitted by a least squares approach? Try both and see which you prefer.

clear;
load('y1_clean.mat'); load('x_clean.mat');
logy=log(y1); logx=log(x);

% Method A: Right line is horizontal, with y-intercept=mean(last 100 points)
% - which is an estimate of the y-asymptote of the curve.
[a0A,a1A,b0A,x0A] = bilinearFitA(logx,logy);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

% Method B: Right line is horizontal, with y-intercept determined by least squares
% fit of all points to the right of the intercept of the left and right lines.
[a0B,a1B,b0B,x0B] = bilinearFitB(logx,logy);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

% plot results
figure;
subplot(211); 
plot(logx,logy,'r.',[min(logx),x0A],a0A+a1A*[min(logx),x0A],'-b',...
    x0A,a0A+a1A*x0A,'gs',[x0A,max(logx)],[b0A,b0A],'-b'); 
grid on; ylabel('log y'); xlabel('log x'); title('Method A')
subplot(212); 
plot(logx,logy,'r.',[min(logx),x0B],a0B+a1B*[min(logx),x0B],'-b',...
    x0B,a0B+a1B*x0B,'gs',[x0B,max(logx)],[b0B,b0B],'-b'); 
grid on; ylabel('log y'); xlabel('log x'); title('Method B')

function [a0,a1,b0,x0] = bilinearFitA(x,y)
%BILINEARFITA Fit two straight lines to a set of x,y data.
% The right line is horizontal, with intercept b0=mean(last 100 y-values)
%   The fitting function is
%   y = a0 + a1*x   if x<=x0
%   y = b0          if x>x0
%   where x0 = (b0-a0)/a1
% x, y should be vectors of equal length.

%% Determine b0 and initial guesses for a0, a1
b0=mean(y(end-99:end));
x1=x(1:floor(length(x)/2));
y1=y(1:floor(length(x)/2));
p = polyfit(x1,y1,1);
a0init=p(2); a1init=p(1);
f0=[a0init;a1init];  % initial guess

% Define function myFun=difference between measured and predicted y
function v=myFun(f)
    %f=[a0;a1];
    a0=f(1); a1=f(2);
    x0=(b0-a0)/a1;
    v=zeros(length(x),1);
    for i=1:length(x)
        if x(i)<=x0
            v(i)=a0+a1*x(i)-y(i);
        else
            v(i)=b0-y(i);
        end
    end
end

%% Find two straight lines
f = lsqnonlin(@myFun,f0);
a0=f(1); 
a1=f(2); 
x0 = (b0-a0)/a1;

end  % end bilinearFitA

function [a0,a1,b0,x0] = bilinearFitB(x,y)
%BILINEARFITB Fit two straight lines to a set of x,y data.
% The right line is horizontal. Its intercept b0 is determined by least
% squares fitting of all the data to the right of x0.
%   The fitting function is
%   y = a0 + a1*x   if x<=x0
%   y = b0          if x>x0
%   where x0 = (b0-a0)/a1
% x, y should be vectors of equal length.

%% Determine initial guesses for a0, a1, b0
x1=x(1:floor(length(x)/2));
y1=y(1:floor(length(x)/2));
x2=x(ceil(length(x)/2):end);
y2=y(ceil(length(x)/2):end);
p = polyfit(x1,y1,1);
a0init=p(2); a1init=p(1);
b0init=mean(y2);
f0=[a0init;a1init;b0init];  % initial guess

% Define function myFun=difference between measured and predicted y
function v=myFun(f)
    %f=[a0;a1;b0];
    a0=f(1); a1=f(2); b0=f(3);
    x0=(b0-a0)/a1;
    v=zeros(length(x),1);
    for i=1:length(x)
        if x(i)<=x0
            v(i)=a0+a1*x(i)-y(i);
        else
            v(i)=b0-y(i);
        end
    end
end

%% Find two straight lines
f = lsqnonlin(@myFun,f0);
a0=f(1); 
a1=f(2); 
b0=f(3); 
x0 = (b0-a0)/a1;

end

Based on your previous comments, I predict that you will not really like either of the plots above, because the left lines in both are not tangent to the red points. They are not tangent because they are minimizing the sum squared error.

You may also think (as I did) that the left lines do not look like they are the best fit lines to the data (the red dots) in the range they span. The lines appear to fit the data from logx=2 to 3.5 much better than the data at the far left, from logx=-0.5 to +0.5. This happens because there are many more points between logx=2 and 3.5 than there are between logx=-0.5 to +0.5. Therefore the left fitted line "tries harder" to fit the data between logx=2 and 3.5 than to fit the far left data.

Good luck.

Answer 2

[William Rose]

@T,

[edit: corrected typo in formula for ]

The figure you included shows two striaight lines (black) fitted to a portion of a curve (blue). You have drawn two dots where the straight lines cross the curve. You want to know how to locate those points. In order to do that, you must explain how you chose the regions of the curve to fit with each straight line.

If you fit a curvilinear data set with two straight lines

if

if

then the two lines intersect at

Here is a figure from an article which I co-authored in which we compared four ways of fitting experimental data: one straight line, blinear i.e. two straight lines (like yours), exponential, and log-log. The different models have dfferent numbers of adjustable parameters, so one "penalizes" models with more adjustable parameters. See the paper (attached) for details, including all equations. The figure from that paper which is most relevant to your question is shown below.

It shows the bilinear fit to data from one subject.

[William Rose]

@T,

The script immersion.m does not run, because file DDATA is missing. I tried renaming the mat file which you attached, data_clean.mat, to DDATA.mat. That led to a different error:

"Unable to resolve the name 'DDATA.group'."

This happens because the data in file data_clean.mat (renamed to DDATA.mat) does not include the variable "DDATA.group" which expected by your script. Please make sure immersion.m works with the files you provide.

[William Rose]

@T,

Your script calls function boundedline(), which is a contribution at the FileExchange. Others helping you may not want to download from the FileExchange, so you may want to avoid calling it. I don't think you need boundedline() for your question. You can just use plot(), or loglog().

[William Rose]

@T,

I uploaded function bilinearFit.m to the File Exchange, here. It fits two straight lines to a set of x,y data. You drew straight lines on a log-log plot in the figure you uploaded. Therefore, if you use bilinearFit(), you should pass the (log x, log y) data to bilinearFit().

[T]

Many thanks @William Rose! I tried the bilinearFit. I also used the test file you have made and updated with my input data. I see this now, the first part of the fit is not exactly a tangent to the curve and hence does not specify the correct slope of the first part of the result.

I want to know the slope of each part of the curve. So a correct fit would be something similar to the pink line in the figure. The second part seems to have a much better fitting curve. I've added the input data here. Could you please suggest how can i achieve this?

[William Rose]

@T,

"I want to know the slope of each part of the curve."

Then differentiate it, and plot the result, dy/dx, versus x. Sort the data first, if is not already in ascending-x order. If the plot of dy/dx versus x looks "noisier" than you like, then compute the derivative with a 2- or 3- or more points wide window. In the example below I estimate the derivative with 1-point wide, 10 points, and 100 points-wide window. Plots 2 and 3 below have the exact same data. The only difference is that the y-scale is zoomed on plot 3.

load('y1_clean.mat'); load('x_clean.mat');
dydx1=diff(y1)./diff(x);                                % 1-point-wide estimate of dy/dx
dydx10=(y1(11:end)-y1(1:end-10))./(x(11:end)-x(1:end-10));   % 10-point-wide estimate
dydx100=(y1(101:end)-y1(1:end-100))./(x(101:end)-x(1:end-100)); % 100-point-wide estimate
% plot results
figure;
subplot(311); plot(x,y1,'r.'); grid on; ylabel('y'); xlabel('x')
subplot(312); 
plot(x(2:end),dydx1,'r.',x(6:end-5),dydx10,'b.',x(51:end-50),dydx100,'g.');
grid on; ylabel('dy/dx'); xlabel('x'); 
legend('1 wide','10 wide','100 wide')
subplot(313); 
plot(x(2:end),dydx1,'r.',x(6:end-5),dydx10,'b.',x(51:end-50),dydx100,'g.'); 
grid on; ylabel('dy/dx'); xlabel('x'); 
legend('1 wide','10 wide','100 wide'); ylim([-1e-4 1e-4])

Good luck.

[T]

Thanks @William Rose! I am still trying to make this work. :)

[William Rose]

@T,

Here's another idea, inspired by your earlier comments.

Fit a horizontal line to the right side: intercept b1=mean(last N points). This is a simple and arbitrary way of estimating the horizontal asymptote.

Fit a straight line a0+a1*logx through the data from logx=0 to 2. This is a simple and arbitrary way of getting a line that is more or less tangent to a chunk of the left hand region of the data, which is what you seem to want to do.

The logx coordinate of the intersection is logx0=(b0-a0)/a1.

load('y1_clean.mat'); load('x_clean.mat');
logy=log(y1); logx=log(x);

% Right line is horizontal, with y-intercept=mean(last N points)
N=100;          % number of points used to estimate right horizontal asymptote
logxlim=[0,2];  % values of logx for left line fit

% Fit left and right lines
p = polyfit(logx(logx>=logxlim(1) & logx<=logxlim(2)),logy(logx>=logxlim(1) & logx<=logxlim(2)),1);
a0=p(2); a1=p(1);               % left line
b0=mean(logy(end-N+1:end));     % right line
logx0=(b0-a0)/a1;               % logx coord of intersection

% display results on console
fprintf('a0=%.2f, a1=%.2f, b0=%.2f, logx0=%.2f\n',a0,a1,b0,logx0)

a0=-6.47, a1=0.99, b0=-3.70, logx0=2.80

% plot results
figure;
plot(logx,logy,'r.',[min(logx),logx0],a0+a1*[min(logx),logx0],'-b',...
    logx0,a0+a1*logx0,'gs',[logx0,max(logx)],[b0,b0],'-b'); 
grid on; ylabel('log y'); xlabel('log x'); title('Tangent and Asymptote Fit')

How's that?