How to perform piece-wise linear regression to determine break point?
185 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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!
![Sample.jpg](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233110/Sample.jpeg)
0 Commenti
Risposte (9)
Steven Lord
il 14 Dic 2021
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)
0 Commenti
Kaiguang Zhao
il 2 Apr 2022
Modificato: Kaiguang Zhao
il 28 Dic 2023
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/950574/image.png)
0 Commenti
Image Analyst
il 14 Dic 2021
@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));
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/834000/image.png)
2 Commenti
Angelavtc
il 14 Dic 2021
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
il 15 Dic 2021
@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.
Star Strider
il 7 Ago 2019
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233116/image.png)
7 Commenti
Star Strider
il 14 Dic 2021
@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
il 24 Feb 2022
I tried your code, but it is not powerful with this data
Can you help me to check it.
John D'Errico
il 7 Ago 2019
Modificato: John D'Errico
il 3 Mag 2021
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233126/image.jpeg)
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:
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
9 Commenti
Bruno Luong
il 25 Feb 2022
Obviously you run on the default solver.
Image Analyst
il 25 Feb 2022
@galaxy, see my answer below where I used your actual data:
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.
Image Analyst
il 8 Ago 2019
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.
![00_Screenshot.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233147/00_Screenshot.png)
0 Commenti
Bruno Luong
il 8 Ago 2019
Another tool is my FEX BSFK
load('Sample2.mat')
BSFK(Sample2(:,1),Sample2(:,2),2,[],[],struct('annimation',1));
![BeakPoint.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/233170/BeakPoint.png)
ans =
struct with fields:
form: 'pp'
breaks: [0 100.0000 200]
coefs: [2×2 double]
pieces: 2
order: 2
dim: 1
0 Commenti
Image Analyst
il 25 Feb 2022
@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));
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/906810/image.png)
0 Commenti
Satoshi Okazaki
il 1 Feb 2023
Modificato: Satoshi Okazaki
il 7 Feb 2023
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).
2 Commenti
Satoshi Okazaki
il 7 Feb 2023
>> [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
il 4 Feb 2024
No. Bogartz provides AN algorithm, not THE algorithm. There are many such possible algorithms. I can think of at least a few.
Vedere anche
Categorie
Scopri di più su Descriptive Statistics in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!