function dirichlet_sample_test ( )

%*****************************************************************************80
%
%% DIRICHLET_SAMPLE_TEST tests DIRICHLET_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2016
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIRICHLET_SAMPLE_TEST\n' );
  fprintf ( 1, '  DIRICHLET_SAMPLE samples the Dirichlet distribution;\n' );
  fprintf ( 1, '  DIRICHLET_MEAN computes the Dirichlet mean;\n' );
  fprintf ( 1, '  DIRICHLET_VARIANCE computes the Dirichlet variance.\n' );

  a(1:n) = [ 0.250, 0.500, 1.250 ];

  check = dirichlet_check ( n, a );

  if ( ~ check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIRICHLET_SAMPLE_TEST - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of components N = %6d\n', n );

  r8vec_print ( n, a, '  PDF parameters A:' );

  mean(1:n) = dirichlet_mean ( n, a );

  variance(1:n) = dirichlet_variance ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF mean, variance:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14f  %14f\n', i, mean(i), variance(i) );
  end

  m2(1:n,1:n) = dirichlet_moment2 ( n, a );
  
  r8mat_print ( n, n, m2, '  Second moment matrix:' );

  for i = 1 : nsample
    [ x(1:n,i), seed ] = dirichlet_sample ( n, a, seed );
  end
 
  xmax(1:n) = r8row_max ( n, nsample, x );
  xmin(1:n) = r8row_min ( n, nsample, x );
  mean(1:n) = r8row_mean ( n, nsample, x );
  variance(1:n) = r8row_variance ( n, nsample, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size = %d', nsample )

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Observed Min, Max, Mean, Variance:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %6d  %14g  %14g  %14g  %14g\n', ...
      i, xmin(i), xmax(i), mean(i), variance(i) );
  end

  return
end
