CODEWITHSUNDEEP
matlabcurve-fitting
Bob Dole
Bob Dole

Reputation: 57

MATLAB Fit a line to a histogram

I am just wondering, how would I go about fitting a line to histogram, using the z-counts as weights? An example of this is shown below (although this post just discusses overlaying multiple plots), taken from Scatter plot with density in Matlab).

My initial thought is to make an array consisting of each pixel from the density plot, repeated n times to make a scatter plot (n == the number of counts), then do a linear polyfit. This seems awfully redundant though.

enter image description here

Upvotes: 0

Views: 2497

Answers (2)

chipaudette
chipaudette

Reputation: 1675

If you have the original data, which is a collection of (x,y) points, you simply do a polyfit on the original data:

p = polyfit(x(:),y(:),1);  %linear fit

That will give you a best fit (in the least-squares sense) to the original data, which is what you want.

If you do not have the original data, and you only have the 2D histogram, the approach that you defined (which basically recreates a facsimile of the original data) will give a similar answer as if you did the polyfit on the original data.

Upvotes: 0

chipaudette
chipaudette

Reputation: 1675

The other approach is to do a weighted least squares solution. You need the (x,y) location of each pixel and the number of counts n within each pixel. Then, I think that you'd do the weighted least-squares this way:

%gather your known data...have x,y, and n all in the same order as each other
A = [x(:) ones(length(x),1)];  %here are the x values from your histogram
b = y(:);  %here are the y-values from your histogram
C = diag(n(:));  %counts from each pixel in your 2D histogram

%Define polynomial coefficients as p = [slope; y_offset]; 

%usual least-squares solution...written here for reference
% b = A*p;               %remember, p = [slope; y_offset];
% p = inv(A'*A)*(A'*b);  %remember, p = [slope; y_offset];

%We want to apply a weighting matrix, so incorporate the weighting matrix
% A' * b = A' * C * A * p;  
p = inv(A' * C * A)*(A' * b);  %remember, p = [slope; y_offset];

The biggest uncertainty for me with this solution is whether the C matrix should be made up of n or n.^2, I can never remember. Hopefully, someone can correct me in the comments, if needed.

Upvotes: 2

Related Questions