Problem of robust fitting using the "robustfit" function

I am using the function "robustfit" to fit a plane(3D) but I have a problem: I do three different calls for this function but I have not the same result those are the calls: date: x, y, z (vectors) call-1: p = robustfit([x y],z) normal = [p(2) p(3) -1]/ norm([p(2) p(3) -1]) normal = 0.5448273 0.8371124 -0.0490510 call-2: p = robustfit([y z],x) normal = [-1 p(2) p(3)]/ norm([-1 p(2) p(3)]) normal = 0.5460477 0.8377283 -0.0065613 call-3: p = robustfit([x z],y) normal = [p(2) -1 p(3) ]/ norm([p(2) -1 p(3)]) normal = 0.5451328 0.8374704 -0.0043365 So how can I know which is the correct normal thank you in advance

Answer 1

[Richard Willey]

Hi Massinissa

In your first example, you are fitting Z as a function of X and Y. In the second you are fitting X as a function of Y and Z. In the last, you are predicting Y as a function of X and Z.

Your choice of model depends on precisely what you are trying to explain. You should have some a priori knowledge that explains

1. Why you are creating a model
2. Why you think that the relationship between the variables should be modeled as a plane (using robust fitting)

I'm attaching some code that shows how to use Principal Component Analysis to perform an Orthogonal Regression. This might be what you're trying to do...

clear all
clc

[CleanX, CleanY] = meshgrid(1:10);
CleanX = reshape(CleanX,100,1);
CleanY = reshape(CleanY,100,1);
CleanZ = 3*CleanX + 4*CleanY;

scatter3(CleanX, CleanY, CleanZ, 'filled')
hold on
foo = fit([CleanX, CleanY], CleanZ, 'poly11')
plot(foo)

%%Add noise vectors to all three dimensions

% Note:  Normally, if we were going to show a regression example, we'd
% create Noisy_Z by adding a noise vector to CleanZ but leave CleanX and
% Clean Y as is.  Here, we're using a fitting technique that is designed to
% create a model where there is noise associated with both the dependent
% and the independent variables.

NoisyX = CleanX + randn(100,1);
NoisyY = CleanY + randn(100,1);
NoisyZ = CleanZ + randn(100,1);

% create a scatter plot

figure1 = figure;
scatter3(NoisyX,NoisyY,NoisyZ, '.k')

Fit a plane to the data using OLS

[foo2, GoF] = fit([NoisyX NoisyY],NoisyZ, 'poly11');

% superimpose the plane on the scatter plot
hold on
h1 = plot(foo2)
set( h1, 'FaceColor', 'r' )

%%Calculate the sum of the Squared Errors 

% Calculate the difference between the observed (NoisyX and NoisyY) and 
% the "known" value for CleanZ

resid1 = CleanZ - foo(NoisyX, NoisyY);

% Square residuals
resid1_sqrd = resid1.^2;

% Take the sum of the squared residuals
SSE_OLS = sum(resid1_sqrd)

% Create textbox
annotation(figure1,'textbox',...
    [0.147443514644351 0.802739726027397 0.305903765690377 0.0931506849315069],...
    'String',{['SSE OLS = ' num2str(SSE_OLS)]},...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 1]);

%%Use Principal Component Analysis to perform an Orthogonal Regression

% PCA is based on centering and rotation.
% PCA rotates the data such that dimension with the greatest amount of
% variance is parallel to the X axis.  The dimension with the second
% largest amount of variance will be parallel to the Y axis.  This operation
% defines a plane.  The direction with the third largest variance will be
% parallel to the Z axis.  This dimension defines a set of residuals which
% are at right angles to the XY plane.

[coeff,score,roots] = princomp([NoisyX NoisyY NoisyZ]);

basis = coeff(:,1:2);
normal = coeff(:,3);
pctExplained = roots' ./ sum(roots)

% Translate the output from PCA back to the original coordinate system

[n,p] = size([NoisyX NoisyY NoisyZ]);
meanNoisy = mean([NoisyX NoisyY NoisyZ],1);
Predicted = repmat(meanNoisy,n,1) + score(:,1:2)*coeff(:,1:2)';

% Generate a fit object that represents the output from the Orthogonal Regression

[foo3, Gof2] = fit([Predicted(:,1) Predicted(:,2)], Predicted(:,3), 'poly11')

h2  = plot(foo3)
set( h2, 'FaceColor', 'b' )

%%Calculate the Sum of the Squared Errors for the Orthogonal Regression
% Calculate residuals
resid2 = CleanZ - Predicted(:,3);

% Square residuals
resid2_sqrd = resid2.^2;

% Take the sum of the squared residuals
SSE_TLS = sum(resid2_sqrd)

annotation(figure1,'textbox',...
    [0.147443514644351 0.802739726027397 0.305903765690377 0.0931506849315069],...
    'String',{['SSE OLS = ' num2str(SSE_OLS)], ['SSE TLS = ' num2str(SSE_TLS)]},...
    'FitBoxToText','off',...
    'BackgroundColor',[1 1 1]);

Answer 2

[Massinissa]

thank you for responding :)

I'm using linear fitting cause my data is a cloud of points en 3D coming from a 3D scanner, this cloud represent a plane which I may get its normal. Presently, I'm using Ordinary Least Squares method to fit my plane(in C++) and because I must have the best accuracy, so I'm trying to use bestfit method

my model is: z =a+bx+cy or y= a+bx+cz or x=a+by+cz (all of those 3 models must, in theory, give a same normal of the plane) but it's not the case:(

[Richard Willey]

Hi Massinissa

From the sounds of things, the Principal Component Analysis based technique is the right way to go. If I understand your problem correctly, you want to identify a plane that best describes the data cloud coming from your scanner.

This isn't really a "fitting" task like regression. Rather, you're trying to describe the sources of variance in the model. PCA will work great for this.

regards,

Richard