MATLAB: Problem of robust fitting using the “robustfit” function

fitfittingMATLABplane fittingrobust fit

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

Best Answer

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]);
Related Question