how to fit segmental regession with matlab

Hello,

there are many threaads about this theme, but I couldnt find the one which fits my problem.

I have set of datapoints with X and Y coordinates. From experience I can say, that they are best represented with two straight lines with an Intersection point.

So, there are two linear regression formulas f1(x) = m1 * x + b1 , f2(x) = m2 * x + b2. How can I get the best matching parameter pairs (with minimal squared error) [m1 m2] and [b1 b2] using a matlab fit?

Answer 1

[Rik]

Actually, you are fitting the intersection point as well. The code below is not guaranteed to result in a continuous function, but you can probably adapt the cost function to get that property.

f1=@(x,b) b(1)*x+b(2);
f2=@(x,b) b(1)*x+b(2);
f_composite=@(x,b) ...
    f1(x,b(1:2)).*double(x>b(5)) + ...
    f2(x,b(3:4)).*double(x<=b(5));
initial_guess=[1 1 1 1 0];

OLS=@(b,x,y,f) sum((f(b,x) - y).^2);
% %adapted Ordinary Least Squares cost function
% OLS_as_char=[...
%     'sum((f(b,x) - y).^2) +'... % OLS part
%     'inf^(any(b<bounds(:,1))-1) +'...%converts true to inf
%     'inf^(any(b>bounds(:,2))-1)'];

opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
% Use 'fminsearch' to minimise the 'OLS' function
fit_val = fminsearch(OLS, initial_guess(:), opts,x,y,f_composite);

[Michael Haitz]

The code looks very promising.

But when I tried to reproduce the above example with x values

x=[0,40,80,100,120,150,170,200];

and y values

y=[1865,2855,3608,4057,4343,4389,4415,4478];

I got the following error messages:

Error using -

Matrix dimensions must agree.

Error in @(b,x,y,f)sum((f(b,x)-y).^2)

Error in fminsearch (line 189)

fv(:,1) = funfcn(x,varargin{:});

Also in the function OLS_as_char 'bounds' is an undefined variable.

Could you please clarify this?

[Rik]

The OLS_as_char is an example of how you can implement bounds and still use this method. It turns out I switched the order of variables in one of the definitions, so it interpreted your data as the parameters.

The code below has this error fixed and should work. I would advise you to also have a look at Bruno's answer.

x=[0,40,80,100,120,150,170,200];
y=[1865,2855,3608,4057,4343,4389,4415,4478];

f1=@(x,b) b(1)*x+b(2);
f2=@(x,b) b(1)*x+b(2);
f_composite=@(x,b) ...
    f1(x,b(1:2)).*double(x>b(5)) + ...
    f2(x,b(3:4)).*double(x<=b(5));
initial_guess=[1 1 1 1 0];

OLS=@(b,x,y,f) sum((f(x,b) - y).^2);
% %adapted Ordinary Least Squares cost function
% OLS_as_char=[...
%     'sum((f(b,x) - y).^2) +'... % OLS part
%     'inf^(any(b<bounds(:,1))-1) +'...%converts true to inf
%     'inf^(any(b>bounds(:,2))-1)'];

opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
% Use 'fminsearch' to minimise the 'OLS' function
fit_val = fminsearch(OLS, initial_guess(:), opts,x,y,f_composite);

clc
fprintf('f1=%.1f*x+%.1f\nf2=%.1f*x+%.1f\n(f=f1 for x>%.1f, otherwise f=f2)\n',fit_val)

Which returns:

f1=13.0*x+2351.8
f2=-612.4*x+-1649.9
(f=f1 for x>-0.1, otherwise f=f2)

[Bruno Luong]

This result is wrong, the break point should be around x=112.

Probably fminsearch do poor job if the initial guess is not well chosen.

[Rik]

It seems a quite poorly conditioned problem for fminsearch to solve. It probably requires way more data to work well with a semi-random initial guesses. With an almost correct initial the result can be quite good:

x=[0,40,80,100,120,150,170,200];
y=[1865,2855,3608,4057,4343,4389,4415,4478];

f1=@(x,b) b(1)*x+b(2);
f2=@(x,b) b(1)*x+b(2);
f_composite=@(x,b) ...
    f1(x,b(1:2)).*double(x>b(5)) + ...
    f2(x,b(3:4)).*double(x<=b(5));
initial_guess=[2 4000 20 2000 100];

OLS=@(b,x,y,f) sum((f(x,b) - y).^2);
% %adapted Ordinary Least Squares cost function
% OLS_as_char=[...
%     'sum((f(b,x) - y).^2) +'... % OLS part
%     'inf^(any(b<bounds(:,1))-1) +'...%converts true to inf
%     'inf^(any(b>bounds(:,2))-1)'];

opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
% Use 'fminsearch' to minimise the 'OLS' function
fit_val = fminsearch(OLS, initial_guess(:), opts,x,y,f_composite);

clc
fprintf('f1=%.1f*x+%.1f\nf2=%.1f*x+%.1f\n(f=f1 for x>%.1f, otherwise f=f2)\n',fit_val)

x_trend=linspace(0,200,100);
figure(1),clf(1)
plot(x,y,'*')
hold on
plot(x_trend,f_composite(x_trend,fit_val))

[Bruno Luong]

fminsearch is not a good algoritm. The few times I have tried it, it won't give satisfied result. But at least it's free.

[Rik]

Maybe another cost function would help, but I suspect it only works well enough if you already have an acceptable approximate.

I think a lot of the functionality in the curve fitting toolbox should have been in the base Matlab, which is why I bought it (and the image processing toolbox).

Answer 2

[Bruno Luong]

Look like this recent thread also answers your question