function walker_sampler_test ( )

%*****************************************************************************80
%
%% walker_sampler_test() tests walker_sampler().
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    11 July 2021
%
%  Author:
%
%    Original C version by Warren Smith.
%    MATLAB version by John Burkardt.
%
  n = 10;
  p = 2.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'walker_sampler_test():\n' );
  fprintf ( 1, '  walker_sampler() creates Walker sample vectors Y and A\n' );
  fprintf ( 1, '  for efficiently sampling a discrete probability vector.\n' );
  fprintf ( 1, '  Test the Walker sampler with a Zipf-type probability.\n' );
%
%  Generate a standard Zipf probability vector for cases 1,...,N,
%  with parameter P.
%
  x = zipf_probability ( n, p );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Zipf probabilities\n' );
  fprintf ( 1, '  for N = %d\n', n );
  fprintf ( 1, '  and parameter P = %g\n', p );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     X(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %14.6g\n', i, x(i+1) );
  end
%
%  For better testing, randomly scramble the probabilities.
%
  i = randperm ( n );
  x(2:n+1) = x(i+1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Randomly permuted X:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     X(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %14.6g\n', i, x(i+1) );
  end
%
%  Build the Walker sampler.
%
  [ y, a ] = walker_build ( n, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Built the sampler\n' );
  fprintf ( 1, '  i Y(i) A(i):\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %4d  %14.6g  %4d\n', i, y(i+1), a(i+1) );
  end
%
%  Prepare to count the frequency of each outcome.
%
  count = zeros ( n + 2, 1 );
%
%  Call the sampler many times.
%
  for i = 1 : 100000
    j = walker_sampler ( n, y, a );
    count(j+1) = count(j+1) + 1;
  end
%
%  Compare normalized sample frequencies to the original probabilities in X.
%
  s = 0.0;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  100000 samples:\n' );
  fprintf ( 1, '  prob   #samples:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %14.6g  %6d\n', x(i+1), count(i+1) );
    expval = x(i+1) * 100000;
    t = expval - count(i+1);
    s = s + t * t / expval;
  end
  s = s / n;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  sumvar = %g (should be about 1)\n', s );

  return
end

