walk_move: Update an ensemble of 'walkers' using the 'walk move'

In svdataman/tonic: Generate and test MCMC output using RW-MH or Goodman-Weare methods

Description

walk_move updates an ensemble of 'walkers' using the 'walk move'

Usage

walk_move(posterior, theta, S = NULL, chatter = 0, ...)

Arguments

posterior	(function) name of the log(densiy) function to sampling from
theta	(array) nwalkers (rows) * M+2 (columns) the current position of each walker
S	(integer) size of the complementary sample
chatter	(integer) how verbose is the output? (0=quiet, 1=normal, 2=verbose)
...	(anything else) other arguments needed for posterior function

Details

A simple implementation of the 'walk move' for the ensemble MCMC sampler proposed by Goodman & Weare (2010).

Value

An array containing the updated positions (in M-dimensional space) of each of the nwalkers walkers. The array is nwalkers (rows) * M+2 (columns). The last two columns list whether the proposed move was rejected or accepted (0 or 1, respectively) and the current (log) posteriod value.

Notes

At input the array <theta> gives the current position of each walker. There are nwalkers rows, one for each walker. Each row has M+2 columns. The first 1:M columns are the position of each walker, an M-dimensional vector. On output the position of each walker is updated. And the M+1 column is assigned 0 (rejected) or 1 (accepted) depending on whether each walker's position is updated this cycle (accept/reject the proposed update position). The M+2 column contains the log(posterior density) at each walker's current position (after update).

Each walker's position is updated in turn. The position of walker j is updated by randomly selecting a subset of walkers from the ensemble (the complementary sample). From this sample we effectively compute the covariance matrix, and use this to define a multivariate Normal, with the mean as the current position of walker j. We then propose a new position for walker j by drawing from this multivariate Normal.

In practice we randomly select M+1 walkers' positions (where M is the number of variables, or dimensions of the problem), from the set of all walkers excluding walker j. With M+1 positions we form a 'simplex'. (Although, by setting the S parameter one can increase the size of the complementary sample if desired, forming an S-polytope.) We then compute the mean position of the complementary sample x and the displacements of each of its walkers from this mean: delta_k = (x_k - <x>). We then produce M+1 univariate, Normal random numbers z_k and compute W = sum_k z_k * delta_k. So W is a random displacement made from a weighted sum of the displacements (of each walker in the complementary sample from their mean), with random weights. (See Goodman & Weare eqn 11.)

This updated position for walker j is then accepted or rejected with a probability that depends on the ratio of the posterior densities at the current and the proposed new position. If the proposal is rejected, walker j remains at its current position. Once every walker has been updated ( j = 1, 2, ..., nwalkers) we move to cycle i+1 and repeat the update step.

The random numbers (the size of the jump z and the uniform variate u used to randomly choose accept/reject) are computed before the loop over walkers begins. This is to make the code a little more clear and efficient.

See Also

gw_sampler, stretch_move

svdataman/tonic documentation built on Aug. 2, 2019, 3:21 p.m.