加权最小二乘回归中的垂直偏移

机器算法验证 r matlab 回归 最小二乘
2022-04-09 20:05:03

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

在此处输入图像描述

然而,垂直偏移最小二乘拟合对异常值(不应该用于模型估计的点)并不鲁棒。因此,我现在正在考虑使用加权垂直偏移最小二乘回归方法。该方法有两个步骤:

  1. 计算将用于线估计的每个点的权重因子;
  2. 在加权最小二乘回归方案中执行垂直偏移。

暂时,我最大的问题来自第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

下图对应上面的脚本: 在此处输入图像描述.

1个回答

完全修改的答案,请参阅历史

从您的链接中获取公式。它包含大量迭代您的输入点的总和。相乘w

i=1nxii=1nwixii=1nyii=1nwiyii=1nxi2i=1nwixi2i=1nxiyii=1nwixiyii=1nyi2i=1nwiyi2n=i=1n1i=1nwi

请注意,我之前建议对坐标加权,但这会导致二阶项为了模拟表示点的多重性的效果(即重复次具有相同的效果),对于在您的点集上迭代的每个总和,您必须恰好有一个.wwwi=3i3wwsum(x_weighted.*y_weighted)B_down

有了这个解决方案,并使用代数数字的精确算术来避免数字问题,二次方程的两个解决方案之一在您提供的示例数据上给出了很好的结果。鉴于在正确计算的情况下只有左右,因此数字问题不应该是一个严重的问题,这与我之前对不正确权重的经验相反。我仍然不知道通常哪种解决方案是正确的,您是否总是可以选择具有正平方根的解决方案,或者您是否必须检查二阶导数的符号。B22