function dirichlet_mix_sample_test ( )

%*****************************************************************************80
%
%% DIRICHLET_MIX_SAMPLE_TEST tests DIRICHLET_MIX_SAMPLE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2016
%
%  Author:
%
%    John Burkardt
%
  comp_num = 2;
  elem_num = 3;
  sample_num = 1000;

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIRICHLET_MIX_SAMPLE_TEST\n' );
  fprintf ( 1, '  DIRICHLET_MIX_SAMPLE samples the Dirichlet Mix distribution;\n' );
  fprintf ( 1, '  DIRICHLET_MIX_MEAN computes the Dirichlet Mix mean;\n' );

  a(1:3,1:2) = [ 0.250, 1.500; ...
                 0.500, 0.500; ...
                 1.250, 2.000 ];

  comp_weight(1:2) = [ 1.0, 2.0 ];

  check = dirichlet_mix_check ( comp_num, elem_num, a, comp_weight );

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

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of elements ELEM_NUM =   %6d\n', elem_num );
  fprintf ( 1, '  Number of components COMP_NUM = %6d\n', comp_num );

  r8mat_print ( elem_num, comp_num, a, '  PDF parameters A(ELEM,COMP):' );

  r8vec_print ( comp_num, comp_weight, '  Component weights:' );

  mean(1:elem_num) = dirichlet_mix_mean ( comp_num, elem_num, a, comp_weight );

  r8vec_print ( elem_num, mean, '  PDF mean:' );

  for sample_i = 1 : sample_num
    [ x(1:elem_num,sample_i), comp, seed ] = dirichlet_mix_sample ( ...
      comp_num, elem_num, a, comp_weight, seed );
  end

  xmax(1:elem_num) = r8row_max ( elem_num, sample_num, x );

  xmin(1:elem_num) = r8row_min ( elem_num, sample_num, x );

  mean(1:elem_num) = r8row_mean ( elem_num, sample_num, x );

  variance(1:elem_num) = r8row_variance ( elem_num, sample_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size = %6d\n', sample_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Observed Min, Max, Mean, Variance:\n' );
  fprintf ( 1, '\n' );

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

  return
end
