Which goodness of fit measure is minimized in robust fit?

I have generated a surface plot of different goodness of fit measures against the fit parameters using robust fit, and noticed that no matter what is the goodness of fit measure (sse, R^2, RMSE..), solution in the minimum does not correspond to the best fit solution given by robust fit.

Perhaps I am comparing fit solutions using a wrong goodness of fit measure?

[Mathieu NOE]

hello

it would help if you'd share your data / code / results / expectations

[mila_po]

Hi Mathieu,

Thank you for your interest. I performed a plane fit to a 3-dimensional data set. "g" and "b" are fit parameters, and out of many possibilities, here I used "sse" (sum of squared errors) as a measure of a fit goodness.

Here is the best fit given by robust fit:

And here is the surface of sse against fit parameters. Black star is the best fit solution, and red star is the minimum:

Here is the snippet of a code. Since I was not sure if gof.sse accounts for robust weights or not, I also calculated weights myself. This is actually another issue, that self-calculated sse does not match sse given as output of fit:

data=readmatrix(file_path);

xdata=data(:,1);
ydata=data(:,2);
zdata=data(:,3);

% Set up fittype and options.
ft = fittype( 'g*b*x+b*y', 'independent', {'x', 'y'}, 'dependent', 'z' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares');
opts.Display = 'Off';
opts.Robust = 'Bisquare';
opts.StartPoint = [0.624 0.888];
K = 4.685; %tuning constant

% Fit model to data.
[fitresult, gof, output] = fit( [xdata, ydata], zdata, ft, opts );
coeffs=coeffvalues(fitresult);
p=numel(coeffs);
b_bf=coeffs(1); %bf = "best fit"
g_bf=coeffs(2);


% fit does not print robust weights, so we calculate them here
res = output.residuals;
[Q,~] = qr(output.Jacobian,0);
h = min(.9999, sum(Q.*Q,2)); %leverage
res_adj = res./sqrt(1-h); %adjusted residuals
res_adj_sorted = sort(abs(res_adj));
s = median(res_adj_sorted(max(1,p):end)) / 0.6745; %median absolute deviation
u=res_adj./(K*s); %standardized residuals
w = (abs(u)<1).*(1-u.^2).^2;
%sse_bf = sum(w.*res.^2); % is not the same as gof.sse ???
sse_bf = gof.sse; 

% Plot fit with data.
figure(1)
pt=plot( fitresult, [xdata, ydata], zdata) 
title(sprintf('b = %.4f \n g = %.4f \n sse = %.4f',coeffs(1),coeffs(2),sse_bf),...
    'FontWeight','Normal','FontSize',12)

% lets scan the surface of sse against fit parameters

g = 0.3:0.012:1.5; 
b = 0.3:0.012:1.5; 

[X,Y]=meshgrid(g,b);
sse=zeros(size(X)); % b x g

for i=1:length(g)
    
    for j=1:length(b)
        
        % parameters are fixed, we only "fit" in order to get fit statistics        
        opts.Lower = [b(j) g(i)];
        opts.StartPoint =[b(j) g(i)];
        opts.Upper = [b(j) g(i)];        
        
        
        [fitresult, gof, output] = fit( [xdata, ydata], zdata, ft, opts );
        
        residuals= output.residuals;
        [Q,~] = qr(output.Jacobian,0);
        h = min(.9999, sum(Q.*Q,2)); %leverage
        residuals_adj = residuals./sqrt(1-h); %adjusted residuals        
        residuals_adj_sorted = sort(abs(residuals_adj));
        s = median(residuals_adj_sorted(max(1,p):end)) / 0.6745; %median absolute deviation 
        stand_res=residuals_adj./(K*s); %standardized residuals
        w = (abs(stand_res)<1).*(1-stand_res.^2).^2;
        %sse(j,i)=sum(w.*residuals.^2); %is not the same as gof.sse ???
        
        sse(j,i) = gof.sse;
        
        
    end
end

sse_min = min(min(sse));
sse_n = sse/sse_min; % for visualization purposes normalize to min value

[x,y]=find(sse == sse_min);
g_min=g(y);
b_min=b(x);

figure(2)
contourf(X,Y,sse_n,'LevelList',1:0.007:1.1,'LineStyle','none')
colormap(flipud(parula))
c=colorbar;
c.Label.String = 'sse/sse_min';

hold on
plot(g_bf,b_bf,'p','MarkerSize',8,'MarkerEdgeColor','k','MarkerFaceColor','k')
hold on
plot(g_min,b_min,'p','MarkerSize',8,'MarkerEdgeColor','r','MarkerFaceColor','r')

xlabel('g','FontSize',18);
ylabel('b','FontSize',18);

[mila_po]

P.S. just to add that even when using self-calculated sse, minimum and best fit solution do not match.

[Mathieu NOE]

hello

yes this is a bit disturbing...why would the "best" fit be different from the sse min value, unless there is somewhere a "hiden" weighting that we are not aware of ?

maybe this is not too critical in the sense that both planes would be very close (especially when you look at the max distance that some experimental points do have with the best fit plane)

unfortunately I don't have the CFT so I cannot run your code and my possibilities to help you are quite limited in this scenario

maybe , try to compare this result with other surface fitting tools , some are available on FEx , like this one :

Fast Polynomial Curve and Surface Fitting - File Exchange - MATLAB Central (mathworks.com)

[mila_po]

Hello Mathieu,

Actually, this issue is irrespective of the fit model. I had another case scenario, where I fitted data with with a straight line. So I would discard the option that this is the issue of surface fitting tool.

Best,

Milana