MATLAB: Find Diverging point of two arrays

dvergenceintersectionslope

Hi, I want to find the diverging point of two arrays. Please see the figure. In this figure the point may be around 44-48, where the arrays diverges fast. I am not able to do.

Best Answer

Alternatively, if similar behavior is always expected, you could fit a piecewise linear model to each curve and compare the slopes of the second segment to the first to see which is the more pronounced change. If it's known the two curves always have same relative shape then probably that breakpoint for the one is going to be pretty good estimate. I posted code some months ago demonstrating solving the breakpoint altho I don't recall the question subect...

OK, here's the link--search did work one only second try to think of what used as catchword...

Piecewise Regression

ADDENDUM

Seems to work reasonably well...I first plotted

x=1:length(d); x=x(:);  % x as observation number
sx=std(x); mx=mean(x);  
xpr=(x-mx)/sx;          % standardize for better numeric properties later
figure
plotyy(x,d(1:2:end),x,d(2:2:end))  % see shapes together; same scaling

and concluded are similar enough to only deal with one specifically as representative of all.

So,

>> b=nlinfit(xpr,d(:,1),@piecewise,[-1 0 0 5])  % estimate breakpoint
b =
 -0.1447    0.6965   -0.0294    6.5596
>> yh=piecewise(b,xpr);   % evaluate the fitted regression
>> close                  % start new figure
>> plot(d(:,1),'x')       % plot original
>> hold on
>> plot(yh,'r')           % overlay the fit

Looks pretty good estimate of where the breakpoint would/should be...so what is it?

>> xBreak=b(3)*sx+mx   % convert the break coefficient estimate to x-coordinates
xBreak =
 48.6640

The piecewise regression function is

>> type piecewise
function y=piecewise(b,x)
  % Piecewise linear regression
  %   y = a1 + b1*x, x<C
  %   y = a2 + b2*x, x>C
  %     = a1 + C*(b1-b2) + b2*x
  % Input coefficients as vector of [a1 b1 C b2]
    ix=x>b(3);  % location> breakpoint
    y(ix)=[b(1) + b(3)*(b(2)-b(4))+b(4)*x(ix)];
    y(~ix)=[b(1)+b(2)*x(~ix)];
    y=y(:);

Given the characteristic shapes are so similar, the "divergence" location is pretty much the location at which the two pieces of the first curve change characteristics. This gives at least a reasonable estimate of where that breakpoint is.

It appears that the second portion is nearly quadratic; I've nor worked on the algebra to see about actually implementing that; you need another condition to tie the coefficients together as above where solve for a2 in terms of a1 and C by equating the two pieces at the intersection point, ie,

a1+b1*C = a2+b2*C --> a2 = a1 + (b1-b2)*C

One possible refinement might be to use the above for an initial estimate the fit the quadratic from floor(xBreak:end) and then look for where that residual peaks to the left of that point. This would rely on the data being quadratic as nearly as seems to be in the sample data, probably.

I just stuck in the actual value...

>> bq=polyfit(x(48:end)',d(48:end,1),2) % fit quadratic over the last
bq =
 -0.0025    0.6033  -24.5031
>> yhq=polyval(bq,x(48:end)');           % evaluate it
>> hold on
>> plot(x(48:end)',yhq,'r')              % the fitted quadratic
>> plot(x(48:end)',yhq-d(48:end,1),'g')  % add the residual
>>

Related Solutions

MATLAB: How to add two different linear fits to the same graph

function y=piecewise(coef,x)
% return piecewise linear fit given [b1,m1,c,m2] as coefficient array and vector x
% c is breakpoint, b1 is intercept of first section, m1,m2 are two segment slopes
    b1=coef(1); m1=coef(2);  % reduce parentheses clutter....

    c=coef(3);  m2=coef(4);
    ix=x>=c;                  % location > breakpoint

    y(ix)=[b1+c*(m1-m2)+m2*x(ix)];
    y(~ix)=[b1+m1*x(~ix)];
    y=reshape(y,size(x));
  end

Example use; contrived example but works well in general with real data in my experience...

x1=1:5; x2=x1+5; x=[x1 x2];
y=[polyval([1 1],x1) polyval([3 -10],x2)]+rand(size(x));
plot(x,y,'*')
hold on
b=nlinfit(x,y,@piecewise,[1 1 4 3])
b =
    1.1746    1.0519    5.5482    3.1284
xh=linspace(x(1),x(end));
yh=piecewise(b,xh);
plot(xh,yh)

ADDENDUM

[dpb wrote earlier comment -- "I've never sat down and done the algebra to write as anything but two linear segments but it would be possible to do the same thing but join a straight line and a quadratic..."]

Well, it dawned on me it isn't difficult to do if one limits the implementation to the one specific case of the quadratic portion being arbitrarily set to one section or the other. Since your data have the nonlinear portion to the right of the breakpoint, I chose to do that case--

function y=piecelinquad(coef,x)
% return piecewise mixed linear/quadratic fit given [b1,m1,c,m2,m3] as
% coefficient array and vector x over which to evaluate function
% c is breakpoint, m1, b1 are slope, intercept of first section, and
% m1,m2 are second segment quadratic, linear coefficients
% The quadratic section as coded is always for x>c, linear for x<=c
%
% dp bozarth  26Aug2018
    b1=coef(1); m1=coef(2);  % reduce parentheses clutter....
    c=coef(3);  p=coef(4:5);
    ix=x>c;                  % location > breakpoint
    y(ix)=polyval([p b1+c*m1],x(ix)-c);
    y(~ix)=b1+m1*x(~ix);
    y=reshape(y,size(x));
  end

Testing on a couple of artificial datasets such as the one used above, it seems to work well in locating the breakpoint and minimizing overall SSE; I noticed there may be a tendency to place the breakpoint somewhat farther inside the linear section as the overall error can be less with a little overshoot at the breakpoint by the quadratic. All these tests were with a small number of points, however, the more extensive datasets you have could be better behaved.

It raises the question, however, of whether, depending on the actual problem constraints of just using smoothing spline or similar?

MATLAB: Data fitting and slop

That's what is known as "piecewise" or "broken-stick" regression model. To incorporate the breakpoint into the model is often called a "change point" regression. It's relatively simple to do, but is a nonlinear problem.

The change point model starts with the broken-stick model, i.e.

Y=b0 + b1(X) + b2(X-C) + e

where

Y is the response variable,
X is the covariate, and
CP is the change point, i.e. where the break occurs.
e is the error term

To formulate a changepoint model into a functional form that can use to estimate the coefficients b and C, write something like

CP=b(4);
if (X < CP) then
  Y = b(1) + b(2)*X;
else
  Y = b(1) + b(2)*X + b(3)*(X-CP);
end

as the functional and use Matlab lsqnonlin

I left CP in as the changepoint to make the code read a little more easily while Matlab uses the array for all the coefficients.

ADDENDUM

That is, the "real" function would look more like--

if (X < b(4)) then    % b4 is the estimated breakpoint value
  Y = b(1) + b(2)*X;
else
  Y = b(1) + b(2)*X + b(3)*(X-b(4));
end

Best Answer

Related Solutions

MATLAB: How to add two different linear fits to the same graph

MATLAB: Data fitting and slop

Related Question