function gamma_sample_test ( )

%*****************************************************************************80
%
%% gamma_sample_test() tests gamma_mean(), gamma_sample(), gamma_variance().
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  test_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'GAMMA_SAMPLE_TEST():\n' );
  fprintf ( 1, '  GAMMA_MEAN() computes the Gamma mean;\n' );
  fprintf ( 1, '  GAMMA_SAMPLE() samples the Gamma distribution;\n' );
  fprintf ( 1, '  GAMMA_VARIANCE() computes the Gamma variance.\n' );

  a_test(1:2) = [ 1.0, 2.0 ];
  b_test(1:2) = [ 3.0, 0.5 ];
  c_test(1:2) = [ 2.0, 0.5 ];

  for test_i = 1 : test_num

    a = a_test(test_i);
    b = b_test(test_i);
    c = c_test(test_i);

    check = gamma_check ( a, b, c );

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

    xmean = gamma_mean ( a, b, c );
    variance = gamma_variance ( a, b, c );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  TEST NUMBER: %6d\n', test_i );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  PDF parameter A =             %14g\n', a );
    fprintf ( 1, '  PDF parameter B =             %14g\n', b );
    fprintf ( 1, '  PDF parameter C =             %14g\n', c );
    fprintf ( 1, '  PDF mean =                    %14g\n', xmean );
    fprintf ( 1, '  PDF variance =                %14g\n', variance );
  
    for i = 1 : nsample
      x(i) = gamma_sample ( a, b, c );
    end

    xmean = mean ( x );
    variance = r8vec_variance ( nsample, x );
    xmax = max ( x );
    xmin = min ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Sample size =     %6d\n', nsample );
    fprintf ( 1, '  Sample mean =     %14g\n', xmean );
    fprintf ( 1, '  Sample variance = %14g\n', variance );
    fprintf ( 1, '  Sample maximum =  %14g\n', xmax );
    fprintf ( 1, '  Sample minimum =  %14g\n', xmin );

  end

  return
end
