与原生最小二乘拟合方案相比,垂直偏移最小二乘拟合有很多优点。下图说明了两者之间的区别,对于这两种方法的更详细的比较,我们参考这里。

然而,垂直偏移最小二乘拟合对异常值(不应该用于模型估计的点)并不鲁棒。因此,我现在正在考虑使用加权垂直偏移最小二乘回归方法。该方法有两个步骤:
- 计算将用于线估计的每个点的权重因子;
- 在加权最小二乘回归方案中执行垂直偏移。
暂时,我最大的问题来自第2步。假设给定权重因子,如何得到公式来估计线的参数?非常感谢!
编辑:
基于@MvG 的善意建议,我在 MATLAB 中实现了该算法:
function line = estimate_line_ver_weighted(pt_x, pt_y,w);
% pt_x x coordinate
% pt_y y coordinate
% w weighting factor
pt_x = pt_x(:);
pt_y = pt_y(:);
w = w(:);
% step 1: calculate n
n = sum(w(:));
% step 2: calculate weighted coordinates
y_square = pt_y(:).*pt_y(:);
x_square = pt_x(:).*pt_x(:);
x_square_weighted = x_square.*w;
y_square_weighted = y_square.*w;
x_weighted = pt_x.*w;
y_weighted = pt_y.*w;
% step 3: calculate the formula
B_upleft = sum(y_square_weighted)-sum(y_weighted).^2/n;
B_upright = sum(x_square_weighted)-sum(x_weighted).^2/n;
B_down = sum(x_weighted(:))*sum(y_weighted(:))/n-sum(x_weighted.*pt_y);
B = 0.5*(B_upleft-B_upright)/B_down;
% step 4: calculate b
if B<0
b = -B+sqrt(B.^2+1);
else
b = -B-sqrt(B.^2+1);
end
% Step 5: calculate a
a = (sum(y_weighted)-b*sum(x_weighted))/n;
% Step 6: the model is y = a + bx, and now we transform the model to
% a*x + b*y + c = 0;
c_ = a;
a_ = b;
b_ = -1;
line = [a_ b_ c_];
结果与我们预期的一样好,如下脚本所示:
%% Procedure 1: given the data
pt_x = [ 692 692 693 692 693 693 750];
pt_y = [ 919 971 1022 1074 1126 1230 1289];
% Procedure 2: draw the point
close all; figure; plot(pt_x,pt_y,'b*');
% Procedure 3: estimate the line based on the weighted vertical offset
% least square method.
weighting = ones(length(pt_x),1);
weighting(end) = 0.01; % we give the last point a low weighting because obvously it is an outlier
myline = estimate_line_ver_weighted(pt_x,pt_y,weighting);
a = myline(1); b = myline(2); c= myline(3);
% Procedure 4: draw the line
x_range = [min(pt_x):0.1:max(pt_x)];
y_range = [min(pt_y):0.1:max(pt_y)];
if length(x_range)>length(y_range)
x_range_corrspond = -(a*x_range+c)/b;
hold on; plot(x_range,x_range_corrspond,'r');
else
y_range_correspond = -(b*y_range+c)/a;
hold on; plot(y_range_correspond,y_range,'r');
end
下图对应上面的脚本:
.