function asa266_test09 ( )

%*****************************************************************************80
%
%% TEST09 tests DIRICHLET_ESTIMATE, DIRICHLET_MEAN, DIRICHLET_VARIANCE, DIRICHLET_SAMPLE.
%
%  Discussion:
%
%    Data is generated by sampling a distribution with known parameters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2010
%
%  Author:
%
%    John Burkardt
%
  elem_num = 3;
  sample_num = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST09\n' );
  fprintf ( 1, '  For a Dirichlet distribution,\n' );
  fprintf ( 1, '  DIRICHLET_SAMPLE samples;\n' );
  fprintf ( 1, '  DIRICHLET_MEAN finds the means;\n' );
  fprintf ( 1, '  DIRICHLET_VARIANCE finds the variances;\n' );
  fprintf ( 1, '  DIRICHLET_ESTIMATE estimates the parameters.\n' );
%
%  Report.
%
  alpha(1:3) = [ 3.22, 20.38, 21.68 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Distribution parameters are:\n' );
  fprintf ( 1, '\n' );
  for elem_i = 1 : elem_num
    fprintf ( 1, '  %4d  %12f\n', elem_i, alpha(elem_i) );
  end

  mean = dirichlet_mean ( elem_num, alpha );

  variance = dirichlet_variance ( elem_num, alpha );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Distribution means, variances are:\n' );
  fprintf ( 1, '\n' );
  for elem_i = 1 : elem_num
    fprintf ( 1, '  %4d  %12f  %12f\n', ...
      elem_i, mean(elem_i), variance(elem_i) );
  end
%
%  Sample the distribution.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of samples is %d\n', sample_num );

  for sample_i = 1 : sample_num

    [ x, seed ] = dirichlet_sample ( elem_num, alpha, seed );

    for elem_i = 1 : elem_num
      x_sample(sample_i,elem_i) = x(elem_i);
    end

  end
%
%  Print some results.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  First few samples:\n' );
  fprintf ( 1, '\n' );

  for sample_i = 1 : min ( sample_num, 10 )
    fprintf ( 1, '  %4d  %12f  %12f  %2f\n', ...
      sample_i, x_sample(sample_i,1:elem_num) );
  end
%
%  Compute means, variances.
%
  mean = r8col_mean ( sample_num, elem_num, x_sample );

  variance = r8col_variance ( sample_num, elem_num, x_sample );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Observed means, variances are:\n' );
  fprintf ( 1, '\n' );
  for elem_i = 1 : elem_num
    fprintf ( 1, '  %4d  %12f  %12f\n', elem_i, mean(elem_i), variance(elem_i) );
  end
%
%  Destroy the values of ALPHA.
%
  alpha(1:elem_num) = 0.0;
%
%  Try to recover the values of ALPHA.
%
  init = 1;
  alpha = [];

  [ alpha, rlogl, v, g, niter, s, epsilon, ifault ] = dirichlet_estimate ( ...
    elem_num, sample_num, x_sample, sample_num,  init, alpha );

  if ( ifault ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIRICHLET_ESTIMATE - Fatal error!\n' );
    fprintf ( 1, '  Error code IFAULT = %d\n', ifault );
    error ( 'DIRICHLET_ESTIMATE - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index, Estimate, Lower Limit, Upper Limit:\n' );
  fprintf ( 1, '\n' );

  for elem_i = 1 : elem_num

    vari = v((elem_i*(elem_i-1))/2+elem_i);
    aminus = alpha(elem_i) - 1.96 * sqrt ( vari );
    aplus = alpha(elem_i) + 1.96 * sqrt ( vari );
    fprintf ( 1, '  %4d  %12f  %12f  %2f\n', ...
      elem_i, alpha(elem_i), aminus, aplus );

  end

  alpha_sum = sum ( alpha(1:elem_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Alpha sum is %f\n', alpha_sum );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NORMALIZED VALUES:\n' );
  fprintf ( 1, '  Index, Estimate, Lower Limit, Upper Limit:\n' );
  fprintf ( 1, '\n' );

  for elem_i = 1 : elem_num
    vari = v((elem_i*(elem_i-1))/2+elem_i);
    aminus = ( alpha(elem_i) - 1.96 * sqrt ( vari ) ) / alpha_sum;
    aplus = ( alpha(elem_i) + 1.96 * sqrt ( vari ) ) / alpha_sum;
    fprintf ( 1, '  %4d  %12f  %12f  %2f\n', ...
      elem_i, alpha(elem_i)/alpha_sum, aminus, aplus );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Log likelikhood function = %f\n', rlogl );

  return
end
