function zp = sample_candidate ( chain_index, chain_num, cr, cr_index, ...
  cr_num, gen_index, gen_num, jumprate_table, jumpstep, limits, pair_num, ...
  par_num, z )

%*****************************************************************************80
%
%% SAMPLE_CANDIDATE generates candidate parameter samples.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2018
%
%  Author:
%
%    Original FORTRAN90 version by Guannan Zhang.
%    MATLAB version by 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 parameters:
%
%    Input, integer CHAIN_INDEX, the chain index.
%    1 <= CHAIN_INDEX <= CHAIN_NUM.
%
%    Input, integer CHAIN_NUM, the total number of chains.
%    3 <= CHAIN_NUM.
%
%    Input, real CR(CR_NUM), the CR values.
%
%    Input, integer CR_INDEX, the index of the chosen CR value.
%    1 <= CR_INDEX <= CR_NUM.
%
%    Input, integer CR_NUM, the total number of CR values.
%    1 <= CR_NUM.
%
%    Input, integer GEN_INDEX, the current generation.
%    1 <= GEN_INDEX <= GEN_NUM.
%
%    Input, integer GEN_NUM, the total number of generations.
%    2 <= GEN_NUM.
%
%    Input, real JUMPRATE_TABLE(PAR_NUM), the jumprate table.
%
%    Input, integer JUMPSTEP, forces a "long jump" every
%    JUMPSTEP generations.
%
%    Input, real LIMITS(2,PAR_NUM), limits for the parameters.
%
%    Input, integer PAIR_NUM, the number of pairs of 
%    crossover chains.
%    0 <= PAIR_NUM.
%
%    Input, integer PAR_NUM, the total number of parameters.
%    1 <= PAR_NUM.
%
%    Input, real Z(PAR_NUM,CHAIN_NUM,GEN_NUM), the Markov chain 
%    sample data.
%
%  Output parameters:
%
%    Output, real ZP(PAR_NUM), a candidate parameter sample.
%
%  Local parameters:
%
%    Input, integer JUMP_DIM(JUMP_NUM), the dimensions in which
%    a jump is to be made.
%
%    Local, integer JUMP_NUM, the number of dimensions in which
%    a jump will be made.  0 <= JUMP_NUM <= PAR_NUM.
%
%    Local, real JUMPRATE, the jump rate.
%

%
%  Used to calculate E following a uniform distribution on (-B,+B).
%  Because B is currently zero, the noise term is suppressed.
%
  b = 0.0;
%
%  Pick pairs of other chains for crossover.
%
  r = zeros ( 2, pair_num );

  for i = 1 : pair_num

    while ( true )

      r2(1,1) = r8_uniform_01_sample ( );
      r2(2,1) = r8_uniform_01_sample ( );

      pair(1:2,1) = floor ( r2(1:2,1) * chain_num ) + 1;

      if ( pair(1,1) ~= pair(2,1) && ...
           pair(1,1) ~= chain_index && ...
           pair(2,1) ~= chain_index )
        break
      end

    end

    r(1:2,i) = pair(1:2,1);

  end
%
%  Determine the jump rate.
%
  [ jump_dim, jump_num, jumprate ] = jumprate_choose ( cr, cr_index, ...
    cr_num, gen_index, jumprate_table, jumpstep, par_num );
%
%  Calculate E in equation 4 of Vrugt.
%
  noise_e = zeros ( par_num, 1 );

  for i = 1 : par_num
    noise_e(i) = b * ( 2.0 * r8_uniform_01_sample ( ) - 1.0 );
  end
%
%  Get epsilon value from multinormal distribution                      
%
  epsilon = zeros ( par_num, 1 );

  av = 0.0;
  sd = 1.0E-10;
  for i = 1 : par_num
    epsilon(i) = r8_normal_sample ( av, sd );
  end
%
%  Generate the candidate sample ZP based on equation 4 of Vrugt.
%
  diff = diff_compute ( chain_num, gen_index, gen_num, jump_dim, jump_num, ...
    pair_num, par_num, r, z );

  zp(1:par_num,1) = z(1:par_num,chain_index,gen_index-1);

  zp(1:par_num) = zp(1:par_num) ...
    + ( 1.0 + noise_e(1:par_num) ) .* jumprate .* diff(1:par_num) ...
    + epsilon(1:par_num);
%
%  Enforce limits on the sample ZP.
%
  zp = sample_limits ( limits, par_num, zp );

  return
end
