How to fix my linear fit model?

142 visualizzazioni (ultimi 30 giorni)
Jenni
Jenni il 3 Apr 2024 alle 14:46
Modificato: AJ il 9 Apr 2024 alle 19:31
Ran in:
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
Cris LaPierre il 3 Apr 2024 alle 15:00
Ran in:
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

Accedi per commentare.

Accedi per rispondere a questa domanda.

Risposta accettata

Bruno Luong
Bruno Luong il 5 Apr 2024 alle 7:40
Modificato: Bruno Luong il 5 Apr 2024 alle 9:00
Ran in:
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
Mathieu NOE il 5 Apr 2024 alle 8:37
Ran in:
+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
Jenni il 5 Apr 2024 alle 8:44
Thank you @Bruno Luong for your help! This tip was exactly the thing I was looking for!

Accedi per commentare.

Più risposte (6)

Cris LaPierre
Cris LaPierre il 3 Apr 2024 alle 15:13
Modificato: Cris LaPierre il 3 Apr 2024 alle 17:20
Ran in:
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
Mathieu NOE il 3 Apr 2024 alle 15:21
Ran in:
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
Bruno Luong il 3 Apr 2024 alle 16:20
Modificato: Bruno Luong il 3 Apr 2024 alle 19:17
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:
Sam Chak
Sam Chak il 3 Apr 2024 alle 16:59
Modificato: Sam Chak il 3 Apr 2024 alle 17:33
Ran in:
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')
Mathieu NOE
Mathieu NOE il 5 Apr 2024 alle 7:02
Ran in:
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
AJ il 9 Apr 2024 alle 19:27
Modificato: AJ il 9 Apr 2024 alle 19:31
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

Categorie

Scopri di più su Interpolation in Help Center e File Exchange

Tag

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by