function [ gr, gr_conv, gr_count ] = gr_compute ( chain_num, gen_index, ...
  gen_num, gr_count, gr_num, gr_threshold, par_num, z )

%*****************************************************************************80
%
%% gr_compute() computes the Gelman Rubin statistics R used to check convergence.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 January 2018
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jasper Vrugt, CJF ter Braak, CGH Diks, Bruce Robinson, James Hyman, 
%    Dave Higdon,
%    Accelerating Markov Chain Monte Carlo Simulation by Differential 
%    Evolution with Self-Adaptive Randomized Subspace Sampling,
%    International Journal of Nonlinear Sciences and Numerical Simulation,
%    Volume 10, Number 3, March 2009, pages 271-288.
%
%  Input:
%
%    integer CHAIN_NUM, the total number of chains.
%    3 <= CHAIN_NUM.
%
%    integer GEN_INDEX, the index of the current generation.
%    1 <= GEN_INDEX <= GEN_NUM.
%
%    integer GEN_NUM, the total number of generations.
%    2 <= GEN_NUM.
%
%    integer GR_COUNT, counts the number of 
%    generations at which the Gelman-Rubin statistic has been computed.
%
%    integer GR_NUM, the number of times the Gelman-Rubin
%    statistic may be computed.
%
%    real GR_THRESHOLD, the convergence tolerance for the
%    Gelman-Rubin statistic.
%
%    integer PAR_NUM, the total number of parameters.
%    1 <= PAR_NUM.
%
%    real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov chain 
%    sample data.
%
%  Output:
%
%    real GR(PAR_NUM,GR_NUM), the Gelman-Rubin R statistic.
%
%    logical GR_CONV, the Gelman-Rubin convergence flag.
%
%    integer GR_COUNT, counts the number of 
%    generations at which the Gelman-Rubin statistic has been computed.
%
  gr_count = gr_count + 1;
  ind0 = floor ( gen_index / 2 );
  rnd0 = ind0;

  for par_index = 1 : par_num

    for chain_index = 1 : chain_num
      mean_chain(chain_index) = ...
        sum ( z(par_index,chain_index,ind0:gen_index) ) / rnd0;
    end

    mean_all = sum ( mean_chain(1:chain_num) ) / chain_num;

    b_var = rnd0 ...
      * sum ( ( mean_chain(1:chain_num) - mean_all ) .^ 2 ) ...
      / ( chain_num - 1 );

    s_sum = 0.0;
    for chain_index = 1 : chain_num
      s = sum ( ( z(par_index,chain_index,ind0:gen_index) ...
        - mean_chain(chain_index) ) .^ 2 );
      s_sum = s_sum + s;
    end
    s_sum = s_sum / ( rnd0 - 1.0 );

    w_var = s_sum / chain_num;

    var = ( ( rnd0 - 1.0 ) * w_var + b_var ) / rnd0;

    gr(par_index,gr_count) = sqrt ( var / w_var );

  end
%
%  Set the convergence flag.
%
  gr_conv = true;

  for par_index = 1 : par_num
    if ( gr_threshold < gr(par_index,gr_count) )
      gr_conv = false;
      break
    end
  end

  if ( gr_conv )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GR_COMPUTE:\n' );
    fprintf ( 1, '  GR convergence at iteration: %d\n', gen_index );
  end

  return
end
