Fit a piece-wise linear function to data: fitting 2 lines to data and finding their corresponding parameters

I have a data series which I know it consists of two parts: two linear sections, the first line is a zero gradient line (y=cte) for x<=a and then for x>=a it is a line with negative slope (y=cx+d).

What I need from the fitted lines is the constant of the first line, the breaking point and the slope of second line.

I don't want to find the breaking point visually and then consider a range for each part of my data and fit a line to each section with polyfit. Because I have plenty of these data series for different experiments with different experimental parameters and I want a routine for that so that I can compare the results I get and don't waste time on each(Not applicable for me at all).

Is there any way with optimization toolbox like lsqcurvefit(which is for non linear functions, and mine isn't) or lsqlin (which I don't get what it is exactly useful for!) or any other commands to carry out the task?

Thanks

[Mona Mahboob Kanafi]

P.S. here is an example data

xdata = [12.2872 12.9804 13.3858 13.6735 13.8967 14.0790 14.2331 14.3667 14.4844 14.5898 14.6851 14.7721 14.8522 14.9263 14.9953 15.0598 15.1204 15.1776 15.2317 15.2829 15.3317 15.3783 15.4227 15.4653 15.5061 15.5453 15.5830 15.6194 15.6545 15.6884 15.7212 15.7529 15.7837 15.8136 15.8426 15.8707 15.8981 15.9248 15.9508 15.9761 16.0008 16.0249 16.0484 16.0714 16.0939 16.1159 16.1374 16.1584 16.1790 16.1992 16.2190 16.2385 16.2575 16.2762 16.2945 16.3126 16.3303 16.3477 16.3647 16.3816 16.3981 16.4143 16.4303 16.4461 16.4616 16.4769 16.4919 16.5067 16.5213 16.5357 16.5499 16.5639 16.5777 16.5913 16.6047 16.6179 16.6310 16.6439 16.6567 16.6692 16.6817 16.6939 16.7061 16.7180 16.7299 16.7416 16.7531 16.7645 16.7758 16.7870 16.7981 16.8090 16.8198 16.8305 16.8411 16.8516 16.8619 16.8722 16.8823 16.8924 16.9023 16.9122 16.9219 16.9316 16.9412 16.9507 16.9600 16.9693 16.9786 16.9877 16.9967 17.0057 17.0146 17.0234 17.0321 17.0408 17.0494 17.0579 17.0663]

[Mona Mahboob Kanafi]

ydata = [-44.8315 -45.4474 -46.6216 -46.5157 -47.2616 -49.7745 -48.4660 -51.8028 -50.5850 -52.1008 -50.4559 -52.6069 -52.8109 -53.1357 -52.6937 -53.4826 -53.3985 -54.2248 -53.6445 -54.5711 -54.2041 -54.6613 -54.1587 -55.8746 -54.3710 -56.3852 -55.0516 -55.9513 -55.3922 -55.9970 -55.7719 -56.0088 -55.5322 -56.6816 -55.3374 -56.7804 -56.0726 -57.0712 -56.1217 -57.0194 -56.6750 -57.4527 -56.5316 -57.5121 -56.3678 -57.6196 -56.9653 -57.4000 -57.2219 -57.8479 -57.2651 -58.2165 -57.6188 -58.1782 -57.5746 -57.7840 -57.6231 -57.8788 -57.9893 -58.5047 -57.9163 -58.0979 -58.1877 -58.7239 -57.9096 -58.6091 -58.2226 -58.8291 -58.5516 -58.8669 -58.5024 -58.8673 -58.5974 -59.3598 -58.7621 -59.2599 -59.0381 -59.2021 -58.9810 -59.4507 -58.9599 -59.6172 -59.2401 -59.4016 -59.2807 -60.0499 -59.0273 -59.5286 -59.4651 -60.0894 -59.3093 -59.5916 -59.0826 -59.8148 -59.0891 -60.0655 -59.3487 -60.1129 -59.6878 -59.7178 -59.6194 -60.1329 -59.5056 -59.9367 -59.5624 -60.1054 -59.9644 -59.9386 -59.8286 -59.9900 -59.6114 -60.0611 -59.7340 -60.1450 -59.6563 -60.0991 -59.8658 -59.9346 -59.6339]

Answer 1

[Teja Muppirala]

You data doesn't really look like it has a constant trend at the beginning. In any case, you should still be able to solve it with LSQCURVEFIT. Here I am making some fake data first.

% Random data...
xdata = linspace(-2,3,101);
ydata = log(abs(10./(10+1i*10.^xdata))) + 0.5*randn(size(xdata));

plot(xdata,ydata)
F = @(B,xdata) min(B(1),B(2)+B(3)*xdata);   %Form of the equation
IC = [max(ydata) max(ydata) 0]; %Initial guess
B = lsqcurvefit( F,IC,xdata,ydata,[min(ydata) -inf -inf],[max(xdata) inf 0]);
hold all;
plot(xdata,F(B,xdata));
a = (B(1) - B(2)) / B(3)
cte = B(1)
c = B(2)
d = B(3)

You may need a better rule to guess your initial conditions, but this general idea should work.

[Matt J]

F = @(B,xdata) min(B(1),B(2)+B(3)*xdata);

This is a non-differentiable fitting function. That could create problems.

[Teja Muppirala]

I think that as long as the error^2 smoothly varies with the parameters, it's no problem that the fitting function itself is not differentiable. For example:

% Original function has a kink in it
x = 0:0.01:10;
y = min(3,x) + 0.2*randn(size(x));

% Fitting function is not differentiable
F = @(B,x) min(B,x);
B_est = lsqcurvefit(F,0,x,y) % True value is B = 3

% Show the results
subplot(211)
plot(x,y,x,F(B_est,x),'r');
legend({'Original Data' char(F)});
title('Fitted function is not differentiable');
subplot(212)
ezplot(@(B) norm( F(B,x) - y ).^2,[0 10]), title('But error^2 is smooth');
hold on;
plot(B_est, norm( F(B_est,x) - y ).^2,'ro')

[Matt J]

I think that as long as the error^2 smoothly varies with the parameters, it's no problem that the fitting function itself is not differentiable.

I have my doubts. LSQCURVEFIT has 2 algorithm options: Levenberg-Marquardt and Trust-Region-Reflective. Levenberg-Marquardt computes the Jacobian of the fitting function. Therefore, the fitting function itself must be differentiable. The trust-region method relies on Hessian computations and therefore requires at least twice differentiability, which will not be satisified for

F = @(B,x) min(B,x)

I don't have the Optimization Toolbox, so I cannot run your code. I'm sure it could work in some cases, though I think there would be an element of luck to it.

Answer 2

[Matt J]

This FEX file looks applicable

http://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation

Answer 3

[Mona Mahboob Kanafi]

Thank you so much for your answers, I'm now using Teja's code, I didn't exactly get what you did in this part : min(B(1),B(2)+B(3)*xdata)

I wrote the following lines and it's working on some case and not working on some other!! I don't know what's the problem with it.

Also, I have this huge difficulty with the start point. I usually select my start points so that they approximately lie on the curve, but this is not working all the time. Can you give me some hints on how to select best starting points?

plot(x,y)

ftobj = fittype(@(y0, x0, c, x) y0 * (x <= x0) + ((c*(x-x0))+y0).* (x > x0))

cfobj = fit(x,y,ftobj,'startpoint',[-30, 6, -4])

hold on;

plot(x,cfobj(x));

[Matt J]

Can you give me some hints on how to select best starting points?

If you can choose an upper and lower portion of the curve safely away from x0, then you could apply POLYFIT to each portion to get an initial estimate c and y0. You could then solve for the intersection of the two pieces to get an initial estimate of x0.