MATLAB: How to fit segmental regession with matlab

Curve Fitting ToolboxMATLABsegmental regression fit

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?

Best Answer

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);

Related Solutions

MATLAB: How to fit implicit function to data

I'll just throw it into the curve fitting toolbox, although any tool would have sufficed for a first try. You should recognize that other tools, such as lsqnonlin from the optimization toolbox, nlinfit from stats TB, etc., would not require a major change to make them solve this too.

z = ones(11,1);
ft = fittype( @(a,b,c,x,y) abs(x).^a + abs(y).^b +c*x.*y, ...
            'independent', {'x', 'y'},'dependent', 'z' );
mdl = fit(data,z,ft,'startpoint',[1 1 1])
     General model:
     mdl(x,y) = abs(x).^a+abs(y).^b+c*x.*y
     Coefficients (with 95% confidence bounds):
       a =       1.868  (1.805, 1.932)
       b =       2.355  (2.27, 2.441)
       c =       2.112  (1.944, 2.279)
       
mdl(data)
ans =
            1
       1.1598
       1.2921
            1
       1.1184
      0.97943
      0.96971
       1.0215
       1.0299
      0.94193
            1

It seems to have worked reasonably. Note the CFT does not provide for bound constraints on the parameters. So if any of a,b,c were going negative, I would have needed to use another tool for the task, or been more creative in my use of fit.

One nice feature of the CFT is it tells me if the parameters are going near zero, in the case of parameters that are desired to bo NOT less than zero. Here, the confidence bounds provided by fit tell me that it would not be considering zero for any of them.

A problem with a fit like this, is it is not really what you probably want, in the sense that it would be nice if it were some sort of orthogonal fit, minimizing the an orthogonal distance of the points to the curve. As well, the fitting tool makes implicit assumptions about the "noise" structrue in the data, that are surely not valid in this case. Remember that z (the target here) was constant at 1. What noise? Don't got no steenkin' noise. Not in z at least.

So, really this is a variation of an errors in variables problem. And one thing you should know about errors in variables problems, is they can be nasty as hell to solve. As well, the result can have biases in the parameters, that there will be specific points that will have considerably more weight in your final result than others, and noise in your data (thus in x and in y) can do strange things.

Regardless, did it work?

syms X Y
fimplicit(abs(X).^mdl.a + abs(Y).^mdl.b + mdl.c*X.*Y - 1)
hold on
grid on
plot(x,y,'ro')

Could be better, could be worse. As I said though, true orthogonal fits are way more difficult to implement than the simple thing I did here.

MATLAB: How to fit this data

Your model function contains some redundant terms that should be merged. Below you find the code to estimate your parameters and how to use that fit. Note that your model returns inf for most parameter sets for the cases where speed is 0. How to deal with that is up to you.

You should note that fminsearch is fairly sensitive to the initial guess, so if you have better ones you might not even have to remove the speed==0. It is also possible to implement bounds by letting the OLS return inf for fit values outside of an allowed range.

%load data from mat file
Map3D=load('Map3D');Map3D=Map3D.Map3D;
Speed = 0:10:60;
Power = 5:10:205;
Altitude = [0 300 600 700];
[S,P,A]=ndgrid(Speed,Power,Altitude);
%remove all speed==0 since that returns a value of inf when evaluated
L_speed_nonzero=S~=0;
[P,S,A,Map3D]=deal(P(L_speed_nonzero),S(L_speed_nonzero),...
    A(L_speed_nonzero),Map3D(L_speed_nonzero));
b_initial_guess=[-10 0.1 -200 -0.05 50 40];
%objective least squares function (requires scalar/vector input)
OLS=@(b,x,y,z,v) sum((MyFun(b,x,y,z) - v).^2);
opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
% Use fminsearch to minimise the 'OLS' function
b_fitted=fminsearch(OLS, b_initial_guess(:), opts,S(:),P(:),A(:),Map3D(:));
clc
fprintf('b_fitted=[')
fprintf('%.4e ',b_fitted)
fprintf('];\n')
%%
%you can calculate the predicted values like this:
%(copied from before, but you should keep it as a variable)
b_fitted=[1.6348e+02 1.0524e-05 -3.2690e-03 6.0934e-08 -1.8714e-04 5.3567e-01 ];
Speed = 0:10:60;
Power = 5:10:205;
Altitude = [0 300 600 700];
[S,P,A]=ndgrid(Speed,Power,Altitude);
predicted_value=MyFun(b_fitted,S,P,A);
%check quality of fit:
real_value=load('Map3D');real_value=real_value.Map3D;
delta=predicted_value-real_value;
delta(i
sinf(delta))=[];
fprintf(['difference between non-inf fitted values and true values ',...
    'is:\nmean= %.2e\nabs max= %.2e\n\n'],mean(delta),max(abs(delta)))
function val=MyFun(b,Power,Speed,Altitude)
% val = Power.^-b(1) ...
%     +b(2) ...
%     +b(3)*Speed.^2 ...
%     +b(4).*Speed ...
%     +b(5) ...
%     +b(6)*Altitude.^2 ...
%     +b(7)*Altitude ...
%     +b(8);
val = Power.^-b(1) ...
    +b(2)*Speed.^2 ...
    +b(3).*Speed ...
    +b(4)*Altitude.^2 ...
    +b(5)*Altitude ...
    +b(6);
end

Best Answer

Related Solutions

MATLAB: How to fit implicit function to data

MATLAB: How to fit this data

Related Question