fit a histogramm to a gaussian- or vice versa

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

One thing you'll need to do is to scale the histogram and bell curve so they have the same area. See my answer to this question:

More Answers (1)

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);

Asked:

on 24 Feb 2012

Edited:

on 25 Oct 2013

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!