function [ fit, z ] = chain_outliers ( chain_num, gen_index, gen_num, ...
  par_num, fit, z )

%*****************************************************************************80
%
%% chain_outliers() identifies and modifies outlier chains during burn-in.
%
%  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.
%    2 <= GEN_INDEX <= GEN_NUM.
%
%    integer GEN_NUM, the total number of generations.
%    2 <= GEN_NUM.
%
%    integer PAR_NUM, the total number of parameters.
%    1 <= PAR_NUM.
%
%    real FIT(CHAIN_NUM,GEN_NUM), the likelihood of
%    each sample.
%
%    real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov
%    chain sample data.
%
%  Output:
%
%    real FIT(CHAIN_NUM,GEN_NUM), the likelihood of
%    each sample.
%
%    real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov
%    chain sample data.
%
  klo = floor ( gen_index / 2 );
  knum = gen_index + 1 - klo;

  avg = zeros ( chain_num, 1 );
  for j = 1 : chain_num
    avg(j) = sum ( fit(j,klo:gen_index) ) / knum;
  end
%
%  Set BEST to be the index of the chain with maximum average.
%
  [ avg_max, best ] = max ( avg );
%
%  Determine the indices of the chains having averages 1/4 "above" 
%  and "below" the average.
%
  avg_sorted = sort ( avg );

  ind1 = 1 + round ( 0.25 * chain_num );
  ind3 = 1 + round ( 0.75 * chain_num );

  q1 = avg_sorted(ind1);
  q3 = avg_sorted(ind3);
  qr = q3 - q1;
%
%  Identify outlier chains, and replace their later samples
%  with values from the "best" chain.
%
  outlier_num = 0;
  for j = 1 : chain_num
    if ( avg(j) < q1 - 2.0 * qr )
      outlier_num = outlier_num + 1;
      z(1:par_num,j,gen_index) = z(1:par_num,best,gen_index);
      fit(j,klo:gen_index) = fit(best,klo:gen_index);
    end
  end

  if ( 0 < outlier_num )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'CHAIN_OUTLIERS:\n' );
    fprintf ( 1, ...
      '  At iteration %d found %d outlier chains,\n', gen_index, outlier_num );
    fprintf ( 1, '  whose indices appear below, and for which samples\n' );
    fprintf ( 1, '  from the chain with the largest log likelihood function,\n' );
    fprintf ( 1, '  index number %d will be substituted.\n', best );

    for j = 1 : chain_num
      if ( avg(j) < q1 - 2.0 * qr )
        fprintf ( 1, '  %d\n', j );
      end
    end

  end

  return
end
