Linear polyfit results don't look like they fit

I have data from two sources that I am comparing via a density scatter plot. I want to find the best linear fit to these data, so I am using polyfit. However, the polyfit results when plotted on top of the density scatter plot don't align very well so I'm questioning the accuracy of the polyfit tool. As you can see from the example figure, the green polyfit line (y = 0.47616x - 0.74845) has a shallower slope and doesn't seem to fit on top of the plotted data very well. (This was all done in R2023a)

edit: as the gracious, kind, gentle John D suggested, there was data outside of the -25:25 limit that was affecting the fit. I was not aware that data existed. That was my issue. Thanks, John! So helpful!

Here is the code:

load('SAR_data.mat','SARucatvalid','SARvcatvalid')
load('AMSR_data.mat','AMSRucatvalid','AMSRvcatvalid')
% create the grid
x1 = -25:0.01:25; y1=x1; y2=-x1;
x2 = zeros(1,5001); 
xGrid = -25:0.1:25;
yGrid = -25:0.1:25;
% determine the polyfit
Ppoly = polyfit(AMSRvcatvalid,SARvcatvalid,1);
yfit = Ppoly(1)*xGrid+Ppoly(2);
Rcoeff = corrcoef(AMSRvcatvalid,SARvcatvalid);
rstr = ['r = ' num2str(Rcoeff(1,2))];
fitstr = ['y = ' num2str(Ppoly(1)) 'x + ' num2str(Ppoly(2))];
% calculate the histogram bins for the density scatter plot
Nv = hist3([AMSRvcatvalid  SARvcatvalid],'Ctrs',{-25:0.1:25 -25:0.1:25});
% make the figure
f0 = figure;
set(gcf,'units','normalized','position',[0 0.15 0.5118 0.6333])
surf(xGrid, yGrid, Nv ); view(2); shading interp
colorbar; hold on; set(gca,'ColorScale','log'); clim([0.1 100])
pp = plot(xGrid, yfit,'g--'); pp.ZData = ones(1,length(pp.XData))*1000;
hold on; p1 = plot(x1,x2,'w--'); hold on; p2 = plot(x2,x1,'w--');
p1.ZData = ones(1,length(p1.XData))*1000; p2.ZData = ones(1,length(p2.XData))*1000;
xlabel('AMSR v'); ylabel('SAR v'); title('Ice Motion v for 8 days in January 2025')
a2 = gca; a2.YDir = 'normal';
text(-20,20,rstr,'color','w','fontsize',15)
text(-20,17,fitstr,'color','g','fontsize',15)
set(gca,'fontsize',15)

[Walter Roberson]

polyfit() for degree 1 is very simple mathematically; I would be quite surprised if the result were "wrong" by more than floating point round-off considerations.

[Alexa Ross]

I agree. I would normally trust the math of the function, but I can't wrap my head around why the linear fit doesn't sit on top of the data. It's not like there are a lot of points on the y=0 line to "pull" the fit to a lower slope.

[John D'Errico]

@Alexa Ross

This is NOT an error in polyfit. It is an error in how you are using it, in the data you gave it, and your assumptions about what is the meaning of a linear regression fit.

However, it seems the emporer does not like being told they are strutting around buck naked. So find someone else to tell you how to not fit your data. I'm willing to tell you the truth, and I did explain the problem and what you were doing wrong, but apparently you don't like hearing the truth, so feel free to make up your own version of the truth.

Answer 1

[Matt J]

Because both your x and y data have random errors, polyfit() is not going to be the best tool. I would recommend instead linear2dFit() from this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-to-fit-plot-conics-and-quadrics