How to fix my linear fit model?
9 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
1 Commento
Cris LaPierre
il 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
Risposta accettata
Bruno Luong
il 5 Apr 2024
Modificato: Bruno Luong
il 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)
2 Commenti
Mathieu NOE
il 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
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)
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)
Più risposte (6)
Cris LaPierre
il 3 Apr 2024
Modificato: Cris LaPierre
il 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')
fit2 = fit(laser(~ind)',foto(~ind)','poly1')
% 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
0 Commenti
Mathieu NOE
il 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
0 Commenti
Bruno Luong
il 3 Apr 2024
Modificato: Bruno Luong
il 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1659066/image.png)
1 Commento
Bruno Luong
il 3 Apr 2024
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
Sam Chak
il 3 Apr 2024
Modificato: Sam Chak
il 3 Apr 2024
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1659121/image.png)
%% 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"])
[Sp, gof] = fit(xx', yy', mdl, 'StartPoint', [9, 3, 19, 0.1, -0.2])
%% Plot result
plot(Sp, x, y), grid on, title('Fit data using Softplus function')
0 Commenti
Mathieu NOE
il 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)
% 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
0 Commenti
AJ
il 9 Apr 2024
Modificato: AJ
il 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
0 Commenti
Vedere anche
Categorie
Scopri di più su Interpolation in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!