using 'fit' to add a trend line to series of scatter points

Hi there,

I am trying to add a trendline to a series of scatter points that represent geographic locations and have an estimate of the goodness of fit.

Here is a simplified version of my data- please can someone help me pin point why/ where it is going wrong. This is the output I get and as you can see this cannot be right from looking at the values on my y-axis

f =

Linear model Poly1:

f(x) = p1*x + p2

Coefficients (with 95% confidence bounds):

p1 = 0.1116 (-0.04946, 0.2726)

p2 = 37.37 (32.47, 42.27)

clear all 
close all 

lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000]
lon= [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000]

h=[30.3 30.45 40.72 40.8];
lat_lon_proportions(h)  %function taken from Jonathen Sullivan - available from https://uk.mathworks.com/matlabcentral/fileexchange/32462-correctly-proportion-a-lat-lon-plot
% plot(lat,lon,'k.')
x=lon
y=lat

hold on
f=fit(y,x,'poly1')
plot(f,y,x)

Answer 1

[Chad Greene]

That is the correct first-order least-squares fit. But I think in your plot the '.' marker on the far left side is coincident with the y axis. The result is it looks like the fit is off. Try using 'o' markers instead, as shown below.

Also, as a general practice to help with debugging these kinds of things, I recommend trying not to change variable names unless there's a really compelling reason to do so. Setting x=lon and y=lat just adds confusion. It's particularly confusing in this case because the syntax plot(f,y,x) plots y on the x axis and x on the y axis, but automatically labels the x data as y, and labels the y data as x. I recommend this rewrite:

lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000];
lon = [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000];

f=fit(lat,lon,'poly1')
plot(f,lat,lon,'o')

xlabel latitude
ylabel longitude

Or, if you want to keep longitudes on the x axis, as is common, this:

lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000];
lon = [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000];

f=fit(lon,lat,'poly1')
plot(f,lon,lat,'o')

xlabel longitude
ylabel latitude

[P_L]

Hi Chad the 2nd option you suggested is the one I am wanting as I need Longitude along the x-axis.

My problem is that the intercept calculated using the fit function is 37.37 which when you look at the graph, cannot be correct. The ouput is below - is there a way to correct this and/ or is there a way to measure the goodness of the fit?

Many Thanks,

f = 

     Linear model Poly1:
     f(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =      0.1116  (-0.04946, 0.2726)
       p2 =       37.37  (32.47, 42.27)

[Chad Greene]

It looks like you got that f with lon and lat swapped. If you want longitude to be the "x" values and latitude to be the "y" values, do this:

lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000];
lon = [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000];
f=fit(lon,lat,'poly1')
f = 
     Linear model Poly1:
     f(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =       1.318  (-0.5842, 3.219)
       p2 =       -23.3  (-100.8, 54.22)

Another way is with polyfit:

P = polyfit(lon,lat,1)
P =
    1.3176  -23.3024

[P_L]

Chad, thank you. But that intercept cannot be right?

Besy wishes ,

P

[Chad Greene]

Yes, -23.3 is the correct y intercept if longitude is on the x axis any latitude is on the y axis. Try it out:

lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000];
lon = [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000];

P = polyfit(lon,lat,1); 

plot(lon,lat,'o')
hold on

% Plot the least-squares trend line: 
lon_array = [0 max(lon)]; 
plot(lon_array,polyval(P,lon_array))

The code above calculates the least-squares fit and then uses polyval to evaluate the fit for the longitude values [0 max(lon]. You will see that the trend line crosses y=-23.3 at x=0. But zoom in really close on the cluster of points in the top right corner, and you'll see that the trend line fits the scattered data.