How to perform piece-wise linear regression to determine break point?

Hello,

I've got a data set that shows two distinct linear segments as shown in the figure below. I'm trying to fit the data to the following equations:

y = m1 + n1 * x for x < b

y = m2 + n2 * x for x > b

At the breakpoint b:

m1 + n1*b = m2 + n2 *b

This means that for x > b :

y = m1 + b(n1 - n2) + n2 * x

The purpose of the regression is to determine the break point b using iterative least square regression, but I'm not sure how to do so in matlab. I've attached the sample data.

x=Sample2(:,1);

y=Sample2(:,2);

Thank you in advance for your time and help!

Answer 1

[Steven Lord]

Consider using the ischange function.

% Making some synthetic data
x = 0:200;
y = zeros(size(x));
y(x <= 100) = x(x <= 100)/2;
y(x > 100) = 50 + 2*(x(x > 100)-100);

% Detect the location of the change
c = ischange(y, 'linear', 'SamplePoints', x);

% Plot the original data using + and the change location using o
plot(x, y, '+')
hold on
plot(x(c), y(c), 'ro', 'MarkerSize', 10)

Answer 2

[Kaiguang Zhao]

Many excellent answers have been posted there. Apparently, numerous algorithms are possible for segmented regression and typically, the results are sensitive to the choice of algorithms. For those who may need a Bayesian alternative for time series changepoint detection, one such Matlab implemenation is available here from this FileExchange entry, which is developed and maintained by me. The algorithm is called BEAST. It can be instantly installed by running eval(webread('http://b.link/rbeast',weboptions('cert',''))). Below is a quick example using a simulated time series:

% Quick installation of BEAST to a temporary path on your local drive
eval(webread('http://b.link/rbeast',weboptions('cert','')))  

% A simulated time series from another quesiton asked in this forum
y = [zeros(1,100) 1:100 99:-1:50  50*ones(1,250)] + 10*rand(1,500); 

% Apply beast to y. Here season='none' indicates that y has no periodic/seasonal component
o =  beast(y, 'season','none')    
printbeast(o)
plotbeast(o)

Here is a summary of the number and locations of the breakpoints (i.e., changepoints) detected:

#####################################################################
#                      Trend  Changepoints                          #
#####################################################################
.-------------------------------------------------------------------.
| Ascii plot of probability distribution for number of chgpts (ncp) |
.-------------------------------------------------------------------.
|Pr(ncp = 0 )=0.000|*                                               |
|Pr(ncp = 1 )=0.000|*                                               |
|Pr(ncp = 2 )=0.000|*                                               |
|Pr(ncp = 3 )=0.914|*********************************************** |
|Pr(ncp = 4 )=0.083|*****                                           |
|Pr(ncp = 5 )=0.002|*                                               |
|Pr(ncp = 6 )=0.000|*                                               |
|Pr(ncp = 7 )=0.000|*                                               |
|Pr(ncp = 8 )=0.000|*                                               |
|Pr(ncp = 9 )=0.000|*                                               |
|Pr(ncp = 10)=0.000|*                                               |
.-------------------------------------------------------------------.
|    Summary for number of Trend ChangePoints (tcp)                 |
.-------------------------------------------------------------------.
|ncp_max    = 10   | MaxTrendKnotNum: A parameter you set           |
|ncp_mode   = 3    | Pr(ncp= 3)=0.91: There is a 91.4% probability  |
|                  | that the trend component has  3 changepoint(s).|
|ncp_mean   = 3.09 | Sum{ncp*Pr(ncp)} for ncp = 0,...,10            |
|ncp_pct10  = 3.00 | 10% percentile for number of changepoints      |
|ncp_median = 3.00 | 50% percentile: Median number of changepoints  |
|ncp_pct90  = 3.00 | 90% percentile for number of changepoints      |
.-------------------------------------------------------------------.
| List of probable trend changepoints ranked by probability of      |
| occurrence: Please combine the ncp reported above to determine    |
| which changepoints below are  practically meaningful              |
'-------------------------------------------------------------------'
|tcp#              |time (cp)                  |prob(cpPr)          |
|------------------|---------------------------|--------------------|
|1                 |199.000000                 |1.00000             |
|2                 |252.000000                 |0.92867             |
|3                 |96.000000                  |0.89042             |
|4                 |471.000000                 |0.01800             |
|5                 |413.000000                 |0.00733             |
|6                 |435.000000                 |0.00692             |
|7                 |483.000000                 |0.00679             |
|8                 |448.000000                 |0.00579             |
|9                 |343.000000                 |0.00204             |
|10                |63.000000                  |0.00154             |
.-------------------------------------------------------------------.

Below is the plot. The trend is fitted using a piecewise polynomial model. Again, as a Bayesian method, BEAST assumes the order of the polynomials for individual segments as uknowns. The orders of the polynomial needed to adequately fit the trend are estimated over time, as depicted iin the tOrder subplot below. The 1st and 4th segments are flat lines, so their estimated polynomial orders are close to zeros.

Answer 3

[Image Analyst]

@Annie Kuang and @Angelavtc below is my code with the sample3 attached data. It also works for your sample2 data.

It finds the place where the left and right lines have the max slope difference. It does not depend on the y value increasing by more than some amount like 2. It automatically finds the kink. And the kink is an actual element that is part of both lines like you wanted in your latest comment.

% Find optimal pair of lines to fit noisy data, one line on left side and one line on right side.  Separation/crossing x value is identified.
clc;	% Clear command window.
clear;	% Delete all variables.
close all;	% Close all figure windows except those created by imtool.
workspace;	% Make sure the workspace panel is showing.
fontSize = 18;
markerSize = 20;
%============================================================================================================
% FIRST CREATE SAMPLE X-Y DATA.  Skip (delete) this part if you already have your x and y data vectors.
% Define number of points.
numPoints = 100;
midIndex = round(numPoints/2);
% Load sample data.
s=load('Sample3.mat')
x = s.Sample3(:, 1);
yNoisy = s.Sample3(:, 2);
% Done defining sample data.
%============================================================================================================
% Plot lines.
subplot(2, 2, 1);
plot(x, yNoisy, 'b.', 'MarkerSize', markerSize);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
title('Noisy Y vs. X data', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
%================================================================================================================================================
% NOW SCAN ACROSS DATA FITTING LEFT AND RIGHT PORTION TO LINES.
% 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,
% so let's use at least 25% of the points to make sure we don't get crazy slopes.
% Initialize structure array
for k = 1 : numPoints
    lineData(k).slopeDifferences = 0;
    lineData(k).line1 = [0,0];
    lineData(k).line2 = [0,0];
end
for k = index1 : index2
    % Get data in left side.
    x1 = x(1:k);
    y1 = yNoisy(1:k);
    % Fit a line through the left side.
    coefficients1 = polyfit(x1, y1, 1); % The slope is coefficients1(1).
    % Get data in right side.
    x2 = x(k+1:end);
    y2 = yNoisy(k+1:end);
    % Fit a line through the left side.
    coefficients2 = polyfit(x2, y2, 1); % The slope is coefficients2(1).
    
    % Compute difference in slopes, and store in structure array along with line equation coefficients.
    lineData(k).slopeDifferences = abs(coefficients1(1) - coefficients2(1));
    lineData(k).line1 = coefficients1;
    lineData(k).line2 = coefficients2;
end
% Find index for which slope difference is greatest.
slopeDifferences = [lineData.slopeDifferences]; % Extract from structure array into double vector of slope differences only
% slope1s = struct2table(lineData.line1); % Extract from structure array into double vector of slopes only
% slope2s = [lineData.line2(1)]; % Extract from structure array into double vector of slopes only
[maxSlopeDiff, indexOfMaxSlopeDiff] = max(slopeDifferences)
% Plot slope differences.
subplot(2, 2, 2);
plot(slopeDifferences, 'b.', 'MarkerSize', markerSize);
xlabel('Index', 'FontSize', fontSize);
ylabel('Slope', 'FontSize', fontSize);
grid on;
caption = sprintf('Slope Differences Maximum at Index = %d, x value = %.2f', indexOfMaxSlopeDiff, x(indexOfMaxSlopeDiff));
title(caption, 'FontSize', fontSize);
% Mark it with a red line.
line([indexOfMaxSlopeDiff, indexOfMaxSlopeDiff], [0, maxSlopeDiff], 'Color', 'r', 'LineWidth', 2);
% Show everything together all on one plot.
% Plot lines.
subplot(2, 2, 3:4);
plot(x, yNoisy, 'b.', 'MarkerSize', markerSize);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
hold on;
% Use the equation of line1 to get fitted/regressed y1 values.
slope1 = lineData(indexOfMaxSlopeDiff).line1(1);
intercept1 = lineData(indexOfMaxSlopeDiff).line1(2);
y1Fitted = slope1 * x + intercept1;
% Plot line 1 over/through data.
plot(x, y1Fitted, 'r-', 'LineWidth', 2);
% Use the equation of line2 to get fitted/regressed y2 values.
slope2 = lineData(indexOfMaxSlopeDiff).line2(1);
intercept2 = lineData(indexOfMaxSlopeDiff).line2(2);
y2Fitted = slope2 * x + intercept2;
% Plot line 2 over/through data.
plot(x, y2Fitted, 'r-', 'LineWidth', 2);
%================================================================================================================================================
% FIND THE CROSSING POINT.  IT IS WHERE THE Y VALUES ARE EQUAL.
% Just set the equations equal to each other and solve for x.  
% So y1Fitted = y2Fitted, which means (slope1 * xc + intercept1) = (slope2 * xc + intercept2).  Solving for xc gives:
xc = (intercept2 - intercept1) / (slope1 - slope2);
y1c = slope1 * xc + intercept1; % Will be the same as y2c.
y2c = slope2 * xc + intercept2; % Will be the same as y1c.
% Mark crossing with a magenta line.
% Draw a line up from the x axis to the crossing point.
line([xc, xc], [0, y1c], 'Color', 'm', 'LineWidth', 2);
% Draw a line over from the y axis to the crossing point.
line([0, xc], [y1c, y1c], 'Color', 'm', 'LineWidth', 1);
caption = sprintf('Data with left and right lines overlaid.  Lines cross at (x,y) = (%.4f, %.4f)', xc, y1c);
title(caption, 'FontSize', fontSize);
message1 = sprintf('Left  Equation: y = %.3f * x + %.3f', slope1, intercept1);
message2 = sprintf('Right Equation: y = %.3f * x + %.3f', slope2, intercept2);
message = sprintf('%s\n%s', message1, message2);
fprintf('%s\n', message);
text(5, 100, message, 'Color', 'r', 'FontSize', 15, 'FontWeight', 'bold', 'VerticalAlignment', 'middle');
uiwait(helpdlg(message));

[Angelavtc]

Dear @Image Analyst, thank you very much for the code and the very detailed explanation. Is there a way to modify it so that the kink is given between a line on the left and a power function (a*x^b+c) on the right? Actually, I am trying to do this, but I am not sure if it is possible with your method since it calculates the differences between slopes.

[Image Analyst]

@Angelavtc, You should probably start your own question for that rather than hijack @Annie Kuang's question. But you can though it would be best I think to specify in advance the index where the change occurs. If you don't know where that is and need to find it, it wouldn't be the two-line fitting algorithm like I used - it would be a lot more complicated.

Answer 4

[Star Strider]

Here, it is easiest to detect the break point first, then do the regression. (This is a simple example of a much more complicated problem. See the Wikipedia article on Segmented regression for an extended discussion.)

Try this:

D = load('Sample2.mat');
x = D.Sample2(:,1);
y = D.Sample2(:,2);

ptdif = diff(y)./diff(x);                       % Point-Wise Slope
L1 = [true; ptdif<2];                           % Logical Vector
L2 = ~L1;                                       % Logical Vector

B1 = polyfit(x(L1), y(L1), 1)
F1 = polyval(B1, x(L1));
B2 = polyfit(x(L2), y(L2), 1)
F2 = polyval(B2, x(L2));

figure
plot(x, y, '.')
hold on
plot(x(L1), F1, '-r')
plot(x(L2), F2, '-g')
hold off
grid

producing:

B1 =
                       0.5      4.65739146117998e-15

B2 =
                         2                      -150

[Annie Kuang]

Thank you so very much for your time and answer.

I was wondering if you might know what to do if the breakpoint is not so obvious (Sample3.mat). While I understand the principle of segmented regressions, I lack the nessessary matlab knowledge to perform such computations on matlab.

The final objective is to obtain a plot like the one below.

Thank you once again for your time and effort! It is greatly appreciated.

[Star Strider]

My pleasure.

Try this:

D = load('Sample3.mat');
x = D.Sample3(:,1);
y = D.Sample3(:,2);
ptdif = gradient(y)./gradient(x);               % Point-Wise Slope
fltr = movmean(ptdif, 10);                      % Moving Average Filter
L1 = [fltr <= 2];                               % Logical Vector
L2 = ~L1;                                       % Logical Vector

B1 = polyfit(x(L1), y(L1), 1)
F1 = polyval(B1, x(L1));
B2 = polyfit(x(L2), y(L2), 1)
F2 = polyval(B2, x(L2));
figure
plot(x, y, '.')
hold on
plot(x(L1), F1, '-r')
plot(x(L2), F2, '-g')
hold off
grid

producing:

B1 =
          1.00104529616725       -0.0696864111498263

B2 =
          4.85583982202447         -163.199110122358

and:

The same essential code works here as well. (I changed the diff call to a gradient call, and added the movmean call.)

[Angelavtc]

Dear @Star Strider could you explain why you considered the number 2 as where the breakpoint would be located? And also why 10 in your movmean? It is not clear to me how to identify those numbers using other data. Thank you very much! :)

[Star Strider]

@Angelavtc — I don’t remember at this point (more than 2½ years later). It likely had to do with the characteristics of the slope of the data.

[Angelavtc]

I see @Star Strider. I think I solved it and thank you very much for answering this old post. However, your Segmented regression makes the lines not necessarily be equal at the breakpoint. Is there a way to correct this? Thank you again, I really appreaciate it.

[Star Strider]

@Angelavtc — I’m not certain that I understand the problem. Again, this is nearly 2½ years old. The points do not have to be continuous at the breakpoint, although if they are not, that might require more involved code.

[galaxy]

Hi @Star Strider

I tried your code, but it is not powerful with this data

Can you help me to check it.

Answer 5

[John D'Errico]

The simple way? Use a tool that can solve the problem. Don't write code to do work when you don't understand the fundamentals behind what you need to do.

Here, I'll use my SLM toolbox to solve the problem for you.

slm = slmengine(x,y,'degree',1,'plot','on','knots',3,'interiorknots','free')
slm = 
  struct with fields:

             form: 'slm'
           degree: 1
            knots: [3×1 double]
             coef: [3×1 double]
            stats: [1×1 struct]
     prescription: [1×1 struct]
                x: [101×1 double]
                y: [101×1 double]
    Extrapolation: 'constant'
>> slm.knots
ans =
            0
          100
          200

I told it to plot the curve and the data, but you do not need to do so. The break point lies in this case, at x==100.

Find SLM for free download here:

https://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling

It uses the optimization toolbox.

Could you have solved the problem using other methods? Well, yes. It is not that difficult. Here is a simple code that does the fit, using nothing more complicated than fminbnd.

function [err,fittedlines] = breakfit(interiorbreak,x,y)
  % objective function to estimate the interior break of a simple broken
  % line fit. fittedlines is a cell array, containing the slope and
  % intercepts of the pair of fitted lines.
  
  % ensure that x and y are column vectors.
  x = x(:);
  y = y(:);
  
  nx = numel(x);
  
  breaks = [min(x),interiorbreak,max(x)];
    
  % which points lie in which interval?
  xbins = discretize(x,breaks);
  
  % write the problem in matrix form
  A = [ones(nx,1),x - breaks(1),(x - breaks(2)).*(xbins == 2)];
  
  % we could use pinv here, but it would be slower then backslash, 
  % and I'll be careful to ensure the problem is not singular.
  coef = A\y;
  
  err = norm(y - A*coef);
  
  % unpack the coefficients so we can convert to a pair of line
  % coefficients. I'll do this in a fairly long form so it might be more
  % comprehensible.
  c1 = coef(1);
  s1 = coef(2);
  s2 = coef(3);
  b1 = breaks(1);
  b2 = breaks(2);
  fittedlines = {[s1,c1 - b1*s1],[s2 + s1,c1 - b1*s1 - b2*s2]};

We use it like this:

dx = max(x) - min(x);
[interiorbreak,fval] = fminbnd(@(b2) breakfit(b2,x,y),min(x) + dx/100,max(x) - dx/100);

interiorbreak
interiorbreak =
          100

So it found the break point as x==100, just as SLM found it. Now, we recover the line coefficients by one more call. Here are the line coefficients it finds for your data (slope, y-intercept).

[~,fittedlines] = breakfit(interiorbreak,x,y)
fittedlines =
  1×2 cell array
    {1×2 double}    {1×2 double}
>> fittedlines{:}
ans =
          0.5   5.6951e-07
ans =
            2         -150

[Knut Inge Hvidsten]

Thanks for your solution. I wonder if your suggestion may contain an error, dropping an "- b1*s1" in the 2nd y-intercept?

fittedlines = {[s1,c1 - b1*s1],[s2 + s1,c1 - b1*s1 - b2*s2]};

[John D'Errico]

It seems you are correct. I've changed my solution to reflect your comment. In the test case, it happened that the error made no impact on the solution, so I never noticed the flaw.

[galaxy]

Hi @John D'Errico

I have no optimization toolbox.

I tried your code, but it is not powerful with this data. Can you help me to check it.

One more problem, how can I get index of interiorbreak value?

Thank you so much

[Bruno Luong]

galaxy here is the result with my FEX

load('C:\Users\Bruno\Desktop\data.mat')
BSFK(x,y,2,[],[],struct('annimation',1))

[galaxy]

@Bruno Luong, thank you so much.

But I have no optimization toolbox, so I can't try your FEX

[Bruno Luong]

BSFK does not require optimization tollbox it support qpas solver that can be found elsewhere, however the author's original link is dead (the link not the author). Alternatively if you tell me your platform and can also provide a proper mex file (pitty TMW does not allow to have mex file in their FEX server).

[galaxy]

@Bruno Luong

Thank you for your feedback.

I tried. But result figure is not same with your result. And there is error.

[Bruno Luong]

Obviously you run on the default solver.

You need to download the qpas solver and add into you matlab path.

[Image Analyst]

@galaxy, see my answer below where I used your actual data:

https://www.mathworks.com/matlabcentral/answers/475306-how-to-perform-piece-wise-linear-regression-to-determine-break-point#answer_904410

You can see it gives you the formulas:

Left Equation: y = 6.226 * x + 0.186

Right Equation: y = -0.237 * x + 59.989

Nothing beyond Base MATLAB is needed for that code.

If you like help I gave you, then you can "Vote" for the answer.

Answer 6

[Image Analyst]

I didn't know there was any existing functions to do it, so (a while ago) I just came up with my own intuitive approach. I just assumed the crossover point would be somewhere in the middle of the data. So I tested every point by fitting a line to the data on the left of it, and fitting a line to the data on the right of it. Then I assumed the "kink" or crossing point would be where the difference in slopes on the right and left sides was the greatest. Attached is my demo. Adapt as desired.

Answer 7

[Bruno Luong]

Another tool is my FEX BSFK

load('Sample2.mat')
BSFK(Sample2(:,1),Sample2(:,2),2,[],[],struct('annimation',1));

ans = 

  struct with fields:

      form: 'pp'
    breaks: [0 100.0000 200]
     coefs: [2×2 double]
    pieces: 2
     order: 2
       dim: 1

Answer 8

[Image Analyst]

@galaxy you don't need anything more than polyfit() and other built-in functions. Look, here I used my piecewise linear fit wit your data. Don't be afraid on the length of the code - most of it is just to make fancy plots to help you understand what's going on. If you want you can take out all the informative plotting stuff and end up with much smaller code. You will see it tells you

Left Equation: y = 6.226 * x + 0.186

Right Equation: y = -0.237 * x + 59.989

% Demo by Image Analyst
% Find optimal pair of lines to fit noisy data, one line on left side and one line on right side.  Separation/crossing x value is identified.
% Initialization Steps.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
markerSize = 4;
%============================================================================================================
% FIRST GET SAMPLE X-Y DATA.
s = load('data.mat')
x = s.x;
yInput = s.y;
% Done defining sample data.
%============================================================================================================
% Define number of points.
numPoints = length(yInput)
midIndex = round(numPoints/2);
% Plot lines.
subplot(2, 2, 1);
plot(x, yInput, 'b.', 'MarkerSize', markerSize);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
title('Noisy Y vs. X data', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Save original y axis range.
yRange = ylim;
%================================================================================================================================================
% NOW SCAN ACROSS DATA FITTING LEFT AND RIGHT PORTION TO LINES.
% 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,
% so let's use at least 25% of the points to make sure we don't get crazy slopes.
% Initialize structure array
for k = 1 : numPoints
    lineData(k).slopeDifferences = 0;
    lineData(k).line1 = [0,0];
    lineData(k).line2 = [0,0];
end
for k = index1 : index2
    % Get data in left side.
    x1 = x(1:k);
    y1 = yInput(1:k);
    % Fit a line through the left side.
    coefficients1 = polyfit(x1, y1, 1); % The slope is coefficients1(1).
    % Get data in right side.
    x2 = x(k+1:end);
    y2 = yInput(k+1:end);
    % Fit a line through the left side.
    coefficients2 = polyfit(x2, y2, 1); % The slope is coefficients2(1).
    
    % Compute difference in slopes, and store in structure array along with line equation coefficients.
    lineData(k).slopeDifferences = abs(coefficients1(1) - coefficients2(1));
    lineData(k).line1 = coefficients1;
    lineData(k).line2 = coefficients2;
end
% Find index for which slope difference is greatest.
slopeDifferences = [lineData.slopeDifferences]; % Extract from structure array into double vector of slope differences only
% slope1s = struct2table(lineData.line1); % Extract from structure array into double vector of slopes only
% slope2s = [lineData.line2(1)]; % Extract from structure array into double vector of slopes only
[maxSlopeDiff, indexOfMaxSlopeDiff] = max(slopeDifferences)
% Plot slope differences.
subplot(2, 2, 2);
plot(slopeDifferences, 'b.', 'MarkerSize', markerSize);
xlabel('Index', 'FontSize', fontSize);
ylabel('Slope', 'FontSize', fontSize);
grid on;
caption = sprintf('Slope Differences Maximum at Index = %d, x value = %.2f', indexOfMaxSlopeDiff, x(indexOfMaxSlopeDiff));
title(caption, 'FontSize', fontSize);
% Mark it with a red line.
line([indexOfMaxSlopeDiff, indexOfMaxSlopeDiff], [0, maxSlopeDiff], 'Color', 'r', 'LineWidth', 2);
% Show everything together all on one plot.
% Plot lines.
subplot(2, 2, 3:4);
plot(x, yInput, 'b.', 'MarkerSize', markerSize);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
hold on;
% Use the equation of line1 to get fitted/regressed y1 values.
slope1 = lineData(indexOfMaxSlopeDiff).line1(1);
intercept1 = lineData(indexOfMaxSlopeDiff).line1(2);
y1Fitted = slope1 * x + intercept1;
% Plot line 1 over/through data.
plot(x, y1Fitted, 'r-', 'LineWidth', 2);
% Use the equation of line2 to get fitted/regressed y2 values.
slope2 = lineData(indexOfMaxSlopeDiff).line2(1);
intercept2 = lineData(indexOfMaxSlopeDiff).line2(2);
y2Fitted = slope2 * x + intercept2;
% Plot line 2 over/through data.
plot(x, y2Fitted, 'r-', 'LineWidth', 2);
ylim(yRange);
%================================================================================================================================================
% FIND THE CROSSING POINT.  IT IS WHERE THE Y VALUES ARE EQUAL.
% Just set the equations equal to each other and solve for x.  
% So y1Fitted = y2Fitted, which means (slope1 * xc + intercept1) = (slope2 * xc + intercept2).  Solving for xc gives:
xc = (intercept2 - intercept1) / (slope1 - slope2);
y1c = slope1 * xc + intercept1; % Will be the same as y2c.
y2c = slope2 * xc + intercept2; % Will be the same as y1c.
% Mark crossing with a magenta line.
% Draw a line up from the x axis to the crossing point.
line([xc, xc], [0, y1c], 'Color', 'm', 'LineWidth', 2);
% Draw a line over from the y axis to the crossing point.
line([0, xc], [y1c, y1c], 'Color', 'm', 'LineWidth', 1);
caption = sprintf('Data with left and right lines overlaid.  Lines cross at (x,y) = (%.4f, %.4f)', xc, y1c);
title(caption, 'FontSize', fontSize);
message1 = sprintf('Left  Equation: y = %.3f * x + %.3f', slope1, intercept1);
message2 = sprintf('Right Equation: y = %.3f * x + %.3f', slope2, intercept2);
message = sprintf('%s\n%s', message1, message2);
fprintf('%s\n', message);
text(20, 35, message, 'Color', 'r', 'FontSize', 15, 'FontWeight', 'bold', 'VerticalAlignment', 'middle');
uiwait(helpdlg(message));

Answer 9

[Satoshi Okazaki]

This is a classic problem. Bogartz 1968 provides the algorithm to determine the breakpoint using iterative least square regression. The implementation is here.

>> x = Sample2(:,1);
>> y = Sample2(:,2);
>> [coef,breakPt,R2] = fitBogartz(x,y,0)

coef =

    0.5000    0.0000    2.0000 -150.0000


breakPt =

  100.0000   50.0000


R2 =

    1.0000

Left equation: y = 0.5 x,

Right equation: y = 2 x - 150,

Breakpoint (x,y) = (100,50).

[Satoshi Okazaki]

For @galaxy data,

>> [coef,breakPt,R2] = fitBogartz(x,y,0)

coef =

    6.1881    0.2267   -0.2436   60.6677


breakPt =

    9.3974   58.3789


R2 =

    0.9996
    

Left equation: y = 6.1881 x + 0.2267,

Right equation: y = -0.2436 + 60.6677,

Breakpoint (x,y) = (9.3974, 58.3789).

[John D'Errico]

No. Bogartz provides AN algorithm, not THE algorithm. There are many such possible algorithms. I can think of at least a few.