# How to fix my linear fit model?

18 views (last 30 days)
Jenni on 3 Apr 2024
Edited: AJ on 9 Apr 2024
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?
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 on 3 Apr 2024
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

Bruno Luong on 5 Apr 2024
Edited: Bruno Luong on 5 Apr 2024
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 on 5 Apr 2024
+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
[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 on 5 Apr 2024
Thank you @Bruno Luong for your help! This tip was exactly the thing I was looking for!

Cris LaPierre on 3 Apr 2024
Edited: Cris LaPierre on 3 Apr 2024
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

Mathieu NOE on 3 Apr 2024
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

Bruno Luong on 3 Apr 2024
Edited: Bruno Luong on 3 Apr 2024
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 on 3 Apr 2024

Sam Chak on 3 Apr 2024
Edited: Sam Chak on 3 Apr 2024
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')

Mathieu NOE on 5 Apr 2024
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

AJ on 9 Apr 2024
Edited: AJ on 9 Apr 2024
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