fit a histogramm to a gaussian- or vice versa
Show older comments
Hello folks.
After making a histogramm, i want to compare it with the bell curve. But i have difficulties to manage a fit.
For calculating the bell curve i'm using this function:
function f = gauss_distribution(x, mu, s)
p1 = -.5 * ((x - mu)/s) .^ 2;
p2 = (s * sqrt(2*pi));
f = exp(p1) ./ p2;
What can i do to make them fit? I'm really frustrated right now. Is there a good way?
Greetings
Accepted Answer
More Answers (1)
Image Analyst
on 25 Feb 2012
One way is to use a simple linear least squares fit. It's probably not the best way since you're fitting the log of the histogram counts instead of the counts so it seems to make the amplitude a little less. There may be a non-linear least squares way that's better. But here is the code, for what it's worth (just copy and paste):
% Demo to fit a histogram to a Gaussian
% Uses a linear least squares fit to the log of the counts.
% Fit is only approximate. The amplitude seems to be a bit low.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables.
workspace; % Make sure the workspace panel is showing.
fontSizeTitle = 20;
fontSizeLabel = 13;
% Make up some sample Gaussian distributed data.
numberOfSamples = 5000;
y = 2*randn(numberOfSamples,1) + 10;
[counts bins] = hist(y(:), 100);
% Plot it.
subplot(3,1,1);
bar(bins, counts, 'BarWidth', 1);
ylabel('Counts', 'FontSize', fontSizeLabel);
xlabel('Bin Value', 'FontSize', fontSizeLabel);
title('Random Sample Data', 'FontSize', fontSizeTitle);
grid on;
xLimits = xlim();
yLimits = ylim();
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]); % Maximize figure.
% Do a least squares fit of the histogram to a Gaussian.
% Assume y = A*exp(-(x-mu)^2/sigma^2)
% Take log of both sides
% log(y) = (-1/sigma^2)*x^2 + (2*mu/sigma^2) + (log(A)-mu^2/sigma^2)
% Which is the same as
% lny = a1*x^2 + a2*x + a3
% Now do the least squares fit.
% Don't include and zero bins in the data because log doesn't like 0.
nonZeroBins = counts > 0;
lny = log(counts(nonZeroBins));
coeffs = polyfit(bins(nonZeroBins), lny, 2)
% Find the Gaussian parameters from the least squares parameters.
sigma = sqrt(-1/coeffs(1))
mu = coeffs(2) * sigma^2/2
amplitude = exp(coeffs(3) + mu^2/sigma^2)
% Now we have all 3 Gaussian parameters,
% so reconstruct the data from the Gaussian model.
xModeled = linspace(bins(1), bins(end), 100);
yModeled = amplitude * exp(-(xModeled - mu).^2 / sigma^2);
% Plot it.
subplot(3,1,2);
plot(xModeled, yModeled, 'r-', 'LineWidth', 2);
ylabel('Counts', 'FontSize', fontSizeLabel);
xlabel('Bin Value', 'FontSize', fontSizeLabel);
title('Modeled Data', 'FontSize', fontSizeTitle);
grid on;
% Set up the x axis the same
xlim(xLimits);
ylim(yLimits);
% Now plot both on the same plot
subplot(3,1,3);
bar(bins, counts, 'BarWidth', 1);
hold on;
plot(xModeled, yModeled, 'r-', 'LineWidth', 2);
ylabel('Counts', 'FontSize', fontSizeLabel);
xlabel('Bin Value', 'FontSize', fontSizeLabel);
title('Both Actual and Modeled Data', 'FontSize', fontSizeTitle);
grid on;
% Set up the x axis the same
xlim(xLimits);
ylim(yLimits);
Categories
Find more on Noncentral t Distribution in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!