How to fix my linear fit model?

Hi all!

I am trying to find a proper way to fit two linear models to my data points. I have attached the code below. Hence, I am not sure if my way is the most practical one when it comes to Matlab. My issue is that I am struggling to set one linear from 0 to 18 and then starting another one from the same point until 25.

How could I fix my code?

Thank you for your help!

My code:

close all
clear all

% Data:
laser = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
foto = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

%Effort to find the proper linear models:
rajaarvo = @(I) 0.1*I.*(I>=0 & I<18)+25*I.*(I>=18 & I<=25);
I = 0:0.1:25;

figure
plot(laser,foto,'.')
hold on
plot(I,rajaarvo(I))
xlabel('Current through laser diode (mA)')
ylabel('Current through foto diode (μA)')
grid on

[Cris LaPierre]

You do not appear to be fitting anything here. You are just scaling the values of I by either 0.1 or 25.

I1 = 0:0.1:18;
I2 = 18:0.1:25;

plot(I1,I1*0.1,'-r',I2, I2*25,'-r')
grid on

Answer 1

[Bruno Luong]

Assuming you know where to split the data for left and right lines:

x = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
y = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

left = x <= 18;
Pl = polyfit(x(left),y(left),1);
right = ~left;
Pr = polyfit(x(right),y(right),1);
Q = Pr-Pl;
xbreak = -Q(2)/Q(1); % roots(Q)
xl = [min(x) xbreak];
xr = [xbreak max(x)];

plot(x, y,'or');
hold on
plot(xl,polyval(Pl,xl),'b',xr,polyval(Pr,xr),'b')
xline(xbreak)

[Mathieu NOE]

+1 !

nice , simple and effective !

you can even estimate the break point with the second derivative of y - no need for manual input

x = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
y = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

% find the break point  
absddy = abs(gradient(gradient(y))); % abs of second derivative 
[m,im] = max(absddy);
xbreak_guess = x(im);

left = x <= 0.9*xbreak_guess;  % take some margin (10% of data removed on LHS of xbreak_guess)
Pl = polyfit(x(left),y(left),1);
right = x >= 1.1*xbreak_guess;  % take some margin (10% of data removed on RHS of xbreak_guess)
Pr = polyfit(x(right),y(right),1);
Q = Pr-Pl;
xbreak=-Q(2)/Q(1)

xbreak = 18.6145

xl = [min(x) xbreak];
xr = [xbreak max(x)];
plot(x, y,'or');
hold on
plot(xl,polyval(Pl,xl),'b')
plot(xr,polyval(Pr,xr),'b')
xline(xbreak)

[Jenni]

Thank you @Bruno Luong for your help! This tip was exactly the thing I was looking for!

Answer 2

[Cris LaPierre]

The simplest way in MATLAB is to use fit. The code would be even simpler if you were fitting a single model to your data, instead of 2.

The elbow is because the two models do not intersect at x=18.

% Data:
laser = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
foto = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

% Fit each part of data
ind = laser<18;
fit1 = fit(laser(ind)',foto(ind)','poly1')

fit1 = 
     Linear model Poly1:
     fit1(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =       0.109  (0.09618, 0.1219)
       p2 =     -0.1406  (-0.2693, -0.01193)

fit2 = fit(laser(~ind)',foto(~ind)','poly1')

fit2 = 
     Linear model Poly1:
     fit2(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =       23.25  (20.7, 25.79)
       p2 =      -424.2  (-477.8, -370.5)

% plot results
plot(laser,foto,'o')
hold on
I = 0:0.1:25;
fity = [fit1(I(I<18));fit2(I(I>=18))];
plot(I,fity,'-r')
hold off

Answer 3

[Mathieu NOE]

hello

simply change your equation for the second part with a x shift

rajaarvo = @(I) 0.1*I.*(I>=0 & I<18)+25*I.*(I>=18 & I<=25);

changed into

xc = 18.5;

rajaarvo = @(I) 0.1*I.*(I>=0 & I<xc)+25*(I-xc).*(I>=xc & I<=25);

full code

% Data:
laser = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
foto = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

%Effort to find the proper linear models:
xc = 18.5;
rajaarvo = @(I) 0.1*I.*(I>=0 & I<xc)+25*(I-xc).*(I>=xc & I<=25);
I = 0:0.1:25;
figure
plot(laser,foto,'*-')
hold on
plot(I,rajaarvo(I))
xlabel('Current through laser diode (mA)')
ylabel('Current through foto diode (μA)')
grid on

Answer 4

[Bruno Luong]

This returns the piecewise linear function that is continuous at the break point.

The correct break point is close to 18.6054 and not 18 according to the FEX BSFK

laser = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
foto = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];



pp=BSFK(laser,foto,2,3,[],struct('sigma',1)) % 2: linear spline with 3 knots

% pp result is attached 
xi=linspace(min(laser),max(laser),513);
plot(laser,foto,'.',xi,ppval(pp,xi),'r');

The fit looks very good:

[Bruno Luong]

More methods are proposed in this thread https://www.mathworks.com/matlabcentral/answers/475306-how-to-perform-piece-wise-linear-regression-to-determine-break-point?s_tid=srchtitle

Answer 5

[Sam Chak]

Hi @Jenni

Upon initial inspection, it appears that the data follows the trend of the ReLU function, which is a piecewise linear function. However, if you prefer a continuous smoother curve, I would suggest using a modified Softplus function to fit the data.

%% Data:
x   = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24];
y   = [2e-2 7e-2 13e-2 2e-1 2.8e-1 3.7e-1 4.6e-1 5.7e-1 6.8e-1 7.9e-1 9e-1 1.03 1.18 1.34 1.51 1.72 2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

%% Processing
xx  = linspace(0, 24, 2401);
yy  = interp1(x, y, xx, 'pchip');

%% Softplus function model
mdl = fittype("a*log(1 + exp(b*(x - c))) + d*x + f", dependent="y", independent="x", coefficients=["a", "b", "c", "d", "f"])

mdl = 
     General model:
     mdl(a,b,c,d,f,x) = a*log(1 + exp(b*(x - c))) + d*x + f

[Sp, gof] = fit(xx', yy', mdl, 'StartPoint', [9, 3, 19, 0.1, -0.2])

Sp = 
     General model:
     Sp(x) = a*log(1 + exp(b*(x - c))) + d*x + f
     Coefficients (with 95% confidence bounds):
       a =       8.846  (8.756, 8.936)
       b =       2.869  (2.841, 2.897)
       c =       18.68  (18.68, 18.68)
       d =      0.1094  (0.1072, 0.1115)
       f =     -0.1585  (-0.1804, -0.1367)

gof = struct with fields:
           sse: 131.0003
       rsquare: 1.0000
           dfe: 2396
    adjrsquare: 1.0000
          rmse: 0.2338

%% Plot result
plot(Sp, x, y), grid on, title('Fit data using Softplus function')

Answer 6

[Mathieu NOE]

yet another simple solution, using fminsearch (the no toolbox solution)

% data
laser = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24];
foto = [2*10^-2 7*10^-2 13*10^-2 2*10^-1 2.8*10^-1 3.7*10^-1 4.6*10^-1 ...
    5.7*10^-1 6.8*10^-1 7.9*10^-1 9*10^-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50];

% Fit the data with bilinearFit
[a0,a1,b0,b1,x0]=bilinearFit(laser,foto);

% Display results on console
fprintf('a0, a1, b0, b1=%.3f, %.3f, %.3f, %.3f, x0=%.3f\n',a0,a1,b0,b1,x0)

a0, a1, b0, b1=-0.218, 0.123, -461.843, 24.934, x0=18.606

% Plot results
figure;
xa=[min(laser),x0]; ya=a0+a1*xa;
xb=[x0,max(laser)]; yb=b0+b1*xb;
plot(laser,foto,'k*',xa,ya,'-r',xb,yb,'-b')
hold on
plot(x0,a0+a1*x0,'g+','LineWidth',1.5,'MarkerSize',7)
grid on; xlabel('X'); ylabel('Y'); title('bilinearFit Test')
legend('Data','a_0+a_1x','b_0+b_1x','x_0,y_0','Location','north')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [a0,a1,b0,b1,x0] = bilinearFit(x,y)
%BILINEARFIT Fit two straight lines to a set of x,y data.
%   The fitting function is
%   y = a0 + a1*x   if x<=x0
%   y = b0 + b1*x  if x>x0
%   where x0 = (a0-b0)/(a1-b1)
%  x, y should be vectors of equal length.

% modified M Noe 22/11/2023


if length(x)~=length(y)
    disp('Error in bilinearFit(x,y): x,y must have equal length.')
    return
elseif min(size(x))>1 || min(size(y))>1
    disp('Error in bilinearFit(x,y): x,y must be vectors.')
    return
end

%% Determine initial guesses for a0, a1, b0, b1
% Sort x, y (in case they are not in order), then divide the data into two
% groups with equal number of elements, then fit straight line to each
% group.  These straight lines are the initial guess for the fit.
[xSorted,I] = sort(x);
ySorted = y(I);
 n = floor(length(x)/2);
% n = floor(length(x)/4);
l = numel(xSorted);
 
x1=xSorted(1:n);
y1=ySorted(1:n);
x2=xSorted(l-n:l);
y2=ySorted(l-n:l);

p = polyfit(x1,y1,1);
a0init=p(2); a1init=p(1);
p = polyfit(x2,y2,1);
b0init=p(2); b1init=p(1);
f0=[a0init;a1init;b0init;b1init];  % initial guess

% Define function myFun=difference between measured and predicted y
function v=myFun(f)
    %f=[a0;a1;b0;b1];
    a0=f(1); a1=f(2); b0=f(3); b1=f(4);
    x0=(a0-b0)/(b1-a1);
    v=zeros(length(x),1);
    for i=1:length(x)
        if x(i)<=x0
            v(i)=a0+a1*x(i)-y(i);
        else
            v(i)=b0+b1*x(i)-y(i);
        end
    end
end

%% Find two straight lines
% f = lsqnonlin(@myFun,f0);
f = fminsearch(@(x) norm(myFun(x)),f0);

a0=f(1); 
a1=f(2); 
b0=f(3); 
b1=f(4);
x0 = (a0-b0)/(b1-a1);

end

Answer 7

[AJ]

Here is how I would have solved the problem, going back to linear regression principles.

function explore_two_line_fit
% explore_two_line_fit - Fit two lines to data.
%  - The first line should be LSE fit to the first section of data points, i.e. the points up to an including X0.
%  - The second line should be LSE fit to the remaining data points
%  - The regrssions lines must intersect at X0 (constraint)
laser = [0 1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 ...
    21 22 23 24]; % aka x
foto = [2e-2 7e-2 13e-2 2e-1 2.8e-1 3.7e-1 4.6e-1 ...
    5.7e-1 6.8e-1 7.9e-1 9e-1 1.03 1.18 1.34 1.51 1.72 ...
    2.06 2.72 13.10 35.40 61.30 87.10 112.03 136.50]; % aka y

x = laser;
y = foto;

% Find optimal intersection point
x0_best = inf;
SSE_best = inf;
for x0 = 17:0.1:19
    [m1,m2,b,SSE,X,Y_hat] = solver(x,y,x0);
    if SSE_best > SSE
        SSE_best = SSE;
        x0_best = x0;
    end
    fprintf('x0=%g, SSE=%g\n', x0, SSE)
end

[m1,m2,b,SSE,X,Y_hat] = solver(x,y,x0_best);

% Plot
hold off
plot(x,y,'o','DisplayName','observations');
hold on
plot(x,Y_hat,'.-','DisplayName','predictions');
hold off

end % function explore_two_line_fit

%---------------------------------------------------------------
function [m1,m2,b,SSE,X,Y_hat] = solver(x,y,x0)
% Contruct the linear regression problem:
%   Y = X*beta
% There are three unknown parameters:
%   m1 - The slope of the first line
%   m2 - The slope of the second line
%   b - The intersection point (in the Y direction)
% The two line intersect at (x0,y0). We don't know x0 yet, but we can determine by trial and error.
% These are used in the following model:
%   y = m1 * (x - x0) + b, for x <= x0
%   y = m2 * (x - x0) + b, for x > x0
n = length(x);
assert(length(y) == n)
p = 3; % The number of parameters being solved for
set1 = x <= x0;
set2 = ~set1;
X = zeros(n,p); % Capital X is a matrix of independent variables
X(set1,1) = x(set1) - x0; % Coefficient m1
X(set2,2) = x(set2) - x0; % Coefficient m2
X(:,3) = ones(n,1); % Intercept, b
Y = y(:); % standard linear regression nomenclature, as column vector
beta = X\Y;
[m1,m2,b] = deal(beta(1),beta(2),beta(3));

% Calculate predicted y
% y_hat = zeros(n,1);
% y_hat(set1) = m1 * (x(set1) - x0) + b;
% y_hat(set2) = m2 * (x(set2) - x0) + b;
Y_hat = X*beta;
% sum square error of predictions
e = Y - Y_hat;
SSE = e' * e;

end % function solver