function histogram ( n )

%*****************************************************************************80
%
%% HISTOGRAM tries 3 ways to compute N normal samples, and histograms results.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of samples.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HISTOGRAM:\n' );
  fprintf ( 1, '  Compute N samples of the normal PDF.\n' );
  fprintf ( 1, '  Compare a histogram of the results to the normal PDF itself.\n' );

  nbin = 41;
%
%  Approximate using the sum of M uniforms.
%
  m = 12;
  y1 = rand ( m, n );
  x1 = sum ( y1 ) - m / 2;
  figure ( 1 )
  clf
  hold on
  hist ( x1, nbin );
  f = hist ( x1, nbin );
  scale = f(floor((nbin+1)/2)) * sqrt ( 2.0 * pi );
  grid on
  title ( sprintf ( '%d normal samples using sum of %d uniforms', n, m ), 'Fontsize', 24 )
  x = linspace ( -4, +4, 101 );
  y = scale * exp ( - x.^2 / 2 ) / sqrt ( 2.0 * pi );
  plot ( x, y, 'r', 'Linewidth', 5 );
  hold off
%
%  Box Muller computes pairs.
%
  nhalf = floor ( ( n + 1 ) / 2 );
  y1 = rand ( nhalf, 1 );
  y2 = rand ( nhalf, 1 );
  y3 = sqrt ( - 2.0 * log ( y1 ) ) .* cos ( 2.0 * pi * y2 );
  y4 = sqrt ( - 2.0 * log ( y1 ) ) .* sin ( 2.0 * pi * y2 );
  x2 = reshape ( [ y3; y4 ], 1, 2 * nhalf );
  x2 = x2(1:n);
  figure ( 2 );
  clf
  hold on
  hist ( x2, nbin );
  f = hist ( x2, nbin );
  scale = f(floor((nbin+1)/2)) * sqrt ( 2.0 * pi );
  grid on
  title ( sprintf ( '%d normal samples using Box-Muller', n ), 'Fontsize', 24 )
  x = linspace ( -4, +4, 101 );
  y = scale * exp ( - x.^2 / 2 ) / sqrt ( 2.0 * pi );
  plot ( x, y, 'r', 'Linewidth', 5 );
  hold off
%
%  RANDN works directly with the normal PDF.
%
  x3 = randn ( n, 1 );
  figure ( 3 );
  clf
  hold on
  hist ( x3, nbin );
  f = hist ( x3, nbin );
  scale = f(floor((nbin+1)/2)) * sqrt ( 2.0 * pi );
  grid on
  title ( sprintf ( '%d normal samples using RANDN', n ), 'Fontsize', 24 )
  x = linspace ( -4, +4, 101 );
  y = scale * exp ( - x.^2 / 2 ) / sqrt ( 2.0 * pi );
  plot ( x, y, 'r', 'Linewidth', 5 );
  hold off

  return
end


