demc_sample.m

A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian co

function [saved_states, saved_logprob] = DEMC_sample(logtarget, n, states, blockindex, ...
    current_logprob, temperature_schedule, gamma_schedule, jitter)
% DEMC_SAMPLE Perform Differential Evolution-Markov chain Monte Carlo sampling
% (see References)
%
% Usage
%
% states = DEMC_sample(logtarget, n,  initial_states, blockindex,...
%       initial_logprob, temperature_schedule, gamma_schedule, jitter)
%
% [states, logprob] = DEMC_sample(logtarget, n,  initial_states, blockindex,...  
%       initial_logprob, temperature_schedule, gamma_schedule, jitter)
%
% INPUTS
%
% logtarget - a function accepting a 1-by-d state vector and returning the log
% probability density for that vector. The density may be unnormalized.
% (DEMC may fail if the target density does not have support over R^d.)
%
% n - number of iterates to sample
%
% initial_states - an m-by-d matrix of states where the rows index Markov chains in
% the population and the columns index the random variables defining the
% state.
%
% OPTIONAL INPUTS (may be [], or may be left off the argument list entirely)
%
% blockindex - a cell array of length b, each of the b cells contains the indexes
% of state vector elements belonging to the b'th block. For DE-MC without blocks,
% set blockindex = {1:d}; this is also the default.
%
% initial_logprob - an m-by-1 vector containing the log probability of the
% initial states; if this has been pre-calculated then this can save time,
% especially if this routine is being called in a loop as part of a hybrid
% sampling scheme.
%
% temperature_schedule - schedule for temperatures used in simulated annealing.
% After the schedule is complete, the temperature is set to 1. An empty schedule
% defaults to a constant temperature of 1.
%
% gamma_schedule - schedule for scale factor for vector differences; the default value
% an unchanging schedule of 2.38/sqrt(2d), which is optimal for d-dimensional
% Gaussian targets. For example, a schedule of 9 iterations at gamma = 2.38/sqrt(2.*d)
% followed by one iteration at 0.99 is specified by
% gamma = [2.38/sqrt(2.*d).*ones(1,9) 0.99]
%
% jitter - 1-by-d dimensional vector of scale factors for the jitter; the
% default value is 1e-5.*ones(1,d)
%
% OUTPUT
%
% states - a m-by-d-by-n array of sampled states.
%
% OPTIONAL OUTPUT
%
% logprob - an m-by-n vector containing the log probabilities of the
% final states of each chain.
%
% References
%
% Cajo F.T. Ter Braak, "A Markov Chain Monte Carlo version of the genetic algorithm
% Differential Evolution: easy Bayesian computing for real parameter spaces"
% Stat Comput (2006) 16:239