DEMC_sample(logtarget, n, states, blockindex, ...

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:239249

%% AUTHOR: Corey Yanofsky
%% $DATE: 08-Dec-2007 00:30:16 $
%% $Revision: 1.3 $
%% DEVELOPED: 7.1.0.246 (R14) Service Pack 3
%% FILENAME: block_DEMC_sample.m
%%
%% This work was supported by grants from Genome Quebec and the National 
%% Science and Engineering Research Council of Canada.

% initialize and set defaults
[m d] = size(states);
if nargin < 8 || isempty(jitter), jitter = 1e-5.*ones(1,d);end
if nargin < 7 || isempty(gamma_schedule), gamma_schedule = 2.38./sqrt(2*d);end
if nargin < 6, temperature_schedule = [];end
if(nargin < 5 ||isempty(current_logprob))
    current_logprob = zeros(m, 1);
    for j = 1:m
        current_logprob(j) = logtarget(states(j,:));
    end
end
if nargin < 4 || isempty(blockindex), blockindex = {1:d};end
nblocks = length(blockindex);

% fill out temperature schedule
temperature_schedule = temperature_schedule(:);
temperature_schedule = [temperature_schedule;ones(n - length(temperature_schedule),1)];

% create gamma index
gamidx = mod(0:(n-1),length(gamma_schedule)) + 1;

saved_states = zeros(m,d,n);
saved_logprob = zeros(m,n);

% sample new states n times...
for i = 1:n
    % for each of the m chains...
    for j = 1:m
        % for each block
        for k = 1:nblocks
            [states(j,:) current_logprob(j)] = DEMC_chain_sample(...
                logtarget,  j, states, current_logprob(j), blockindex{k},...
                temperature_schedule(i), gamma_schedule(gamidx(i)), jitter);
        end % blocks
    end % chains
    saved_states(:,:,i) = states;
    saved_logprob(:,i) = current_logprob;
end % iterates

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [new_state new_logprob] = DEMC_chain_sample(logtarget,i,states,current_logprob, ...
    blockindex, temperature, gamma, jitter)
%
% This function does the actual Metropolis sampling.
%
% [new_state, new_logprob] = DEMC_chain_sample(logtarget,i,states,current_logprob, gamma, 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 over R^d.)
%
% 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.
%
% i - the index of the Markov chain being updated
% 
% current_logprob - log probability of the state being updated; if this has
% been pre-calculated (for example, in a previous iteration of DE-MC), then 
% this saves time
%
% blockindex - array containing the indices of the elements of the state
% vectors being updated
%
% temperature - simulated annealing temperature parameter
%
% gamma - scale factor for DE-MC; the default value is 2.38/sqrt(2d), which
% is optimal for d-dimensional Gaussian targets
%
% jitter - 1-by-d dimensional vector of scale factors for the jitter; the
% default value is 1e-5.*ones(1,d)
%
% OUTPUT
%
% new_state - a 1-by-d dimensional state vector sampled by DE-MC.
%
% new_logprob - the log probability of new_state. This will probably not be
% useful when the DE-MC is providing a conditional update as part of a 
% hybrid Markov chain; however, if DE-MC is being used for all 
% updates, then this can be fed back into the next iteration, saving
% computation time.

[m,d] = size(states);

% sample two states other than i uniformly without replacement
n1 = ceil(rand*(m-1));
n2 = ceil(rand*(m-2));
if(n1 >= i) 
    n1 = n1 + 1;
    if(n2 >= i), n2 = n2 + 1;end
    if(n2 >= n1), n2 = n2 + 1;end
else
    if(n2 >= n1), n2 = n2 + 1;end
    if(n2 >= i), n2 = n2 + 1;end
end

% generate proposal via differential evolution rule
difference_vector = zeros(1,d);
difference_vector(blockindex) = gamma.*(states(n1,blockindex) - states(n2,blockindex)) + ...
    jitter(blockindex).*randn(1,length(blockindex));
proposal = states(i,:) + difference_vector;

proposal_logprob = logtarget(proposal);

% Metropolis accept/reject step
if(temperature*log(rand) < proposal_logprob - current_logprob)
    new_state = proposal;
    new_logprob = proposal_logprob;
else
    new_state = states(i,:);
    new_logprob = current_logprob;
end

% Created with NEWFCN.m by Frank Gonzlez-Morphy  
% Contact...: [e-mail]  
% ===== EOF ====== [DEMC_sample.m] ======