Perpendicular offset least square fitting has a lot of advantages compared to the native least square fitting scheme. The following figure illustrates the difference between there, and for a more detailed comparison of these two methods, we refer to here.
Perpendicular offset least square fitting, however, is not robust to outliers( points that are not supposed to be used for model estimation). Therefore, I am now considering to use a weighted perpendicular offset least square regression method. The method has two steps:
- Calculate the weighting factor for each points that are going to be used for line estimation;
- Perform perpendicular offset in a weighted least square regression scheme.
For the time being, my biggest problem comes from step 2. Suppose the weighting factors are given, how can I get the formula to estimate the parameters of the line? Many thanks!
Based on the kind suggestion of @MvG I have implemented the algorithm in 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_];
The result is as good as we can expect, which is illustrated in the following script:
%% 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
The following figure corresponds to the above script: .