How can I fit multiple lines through a common y-intercept?

4 次查看（过去 30 天）
Michael McDonald
Michael McDonald 2020-10-22
编辑： Michael McDonald 2020-10-23
I would like to linearly fit N sets of data to a functional form yi = mi * x + b, such that the lines have individual slopes mi but are required to share a single y-intercept b.
I've included some example code below for a simple case of two vectors y1 and y2 that perfectly lie on lines that would otherwise intercept the y-axis at y=1 and y=2. In this case with identical numbers of points, all perfectly on a line already, I'm pretty sure the "correct" answer would solve to a shared intercept on the black cross at y=1.5. For this simple case I could brute force it by looping through a bunch of y-intercepts, fitting to all of them, and minimizing the residuals, but is there a more elegant linear algebra way to solve this?
Thanks!
%% Initialization
x = [1:5];
m1 = 1;
b1 = 1;
y1 = m1*x + b1;
m2 = 2;
b2 = 2;
y2 = m2*x + b2;
%% Fitting
M1B1 = polyfit(x,y1,1);
M2B2 = polyfit(x,y2,1);
%% Plotting
% figure()
plot(x,y1,'bo',x,y2,'ro',[0 x],[b1 y1],'b--',[0 x],[b2 y2],'r--',0,(b1+b2)/2,'k+')
xlim([-1 5])
ylim([0 10])
title({'What''s the best red and blue line','to fit the data and share a y-intercept?'})

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2020-10-22
编辑：Matt J 2020-10-22
p=[x(:).^[1,0] ,0*x(:); 0*x(:),x(:).^[0,1]] \[y1(:);y2(:)];
m1=p(1)
b=p(2)
m2=p(3)
  2 个评论
Matt J
Matt J 2020-10-22
编辑：Matt J 2020-10-22
For N lines,
Y=vertcat(y1(:),y2(:),...,yN(:));
X=kron(speye(N),x(:));
X(:,N+1)=1;
MB=X\Y;
m=MB(1:N); %vector of N slopes
b=MB(N+1);
Michael McDonald
Michael McDonald 2020-10-23
Thanks Matt! That solves the problem very elegantly indeed. I pasted some code and a picture illustrating the effect, and will mark the question as solved.
%% Initialization
x = [-3:3];
noise = 0.5;
b0 = 1
m1 = 1;
y1 = m1*x + b0 + noise*(rand(size(x))-.5);
m2 = 2;
y2 = m2*x + b0 + noise*(rand(size(x))-.5);
%% Fitting
% Individual fits -- not satisfactory, won't share y-intercept
fit1 = polyfit(x,y1,1)
fitm1 = fit1(1); fitb1 = fit1(2);
fit2 = polyfit(x,y2,1)
fitm2 = fit2(1); fitb2 = fit2(2);
% Matt J's addition: forces a shared y-intercept
p=[x(:).^[1,0] ,0*x(:); 0*x(:),x(:).^[0,1]] \[y1(:);y2(:)];
bestm1=p(1)
b=p(2)
bestm2=p(3)
%% Plotting
figure()
subplot(1,3,1)
plot(x,y1,'bo',x,y2,'ro',x,fitm1*x + fitb1,'b--',x,fitm2*x + fitb2,'r--',0,b0,'k*')
title({'Here are simple individual fits','to these two datasets'})
subplot(1,3,2)
plot(x,y1,'bo',x,y2,'ro',x,fitm1*x + fitb1,'b--',x,fitm2*x + fitb2,'r--',0,b0,'k*')
xlim([-.1 .1])
ylim([b0-.25 b0+.25])
title({'Zoom in: We''d like them','to share a y-intercept'})
subplot(1,3,3)
plot(x,y1,'bo',x,y2,'ro',x,bestm1*x+b,'b-',x,bestm2*x+b,'r-',0,b0,'k*',0,b,'ko')
xlim([-.1 .1])
ylim([b0-.25 b0+.25])
title('Much better!')

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Least Squares 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by