Slice Sampling - Doubling MCMC

class pints.SliceDoublingMCMC(x0, sigma0=None)

Implements Slice Sampling with Doubling, as described in [1].

This is a univariate method, which is applied in a Slice-Sampling-within-Gibbs framework to allow MCMC sampling from multivariate models.

Generates samples by sampling uniformly from the volume underneath the posterior (\(f\)). It does so by introducing an auxiliary variable (\(y\)) and by defining a Markov chain.

If the distribution is univariate, sampling follows:

1. Calculate the pdf (\(f(x0)\)) of the current sample (\(x0\)).
2. Draw a real value (\(y\)) uniformly from (0, f(x0)), defining a horizontal “slice”: \(S = {x: y < f (x)}\). Note that \(x0\) is always within S.
3. Find an interval (\(I = (L, R)\)) around \(x0\) that contains all, or much, of the slice.
4. Draw a new point (\(x1\)) from the part of the slice within this interval.

If the distribution is multivariate, we apply the univariate algorithm to each variable in turn, where the other variables are set at their current values.

This implementation uses the “Doubling” method to estimate the interval \(I = (L, R)\), as described in [1] Fig. 4. pp.715 and consists of the following steps:

1. \(U \sim uniform(0, 1)\)
2. \(L = x_0 - wU\)
3. \(R = L + w\)
4. \(K = p\)

5. while \(K > 0\) and \({y < f(L) or y < f(R)}\):

   1. \(V \sim uniform(0, 1)\)
   2. if \(V < 0.5\), then \(L = L - (R - L)\) else, \(R = R + (R - L)\)

6. \(K = K - 1\)

Intuitively, the interval I is estimated by expanding the initial interval by producing a sequence of intervals, each twice the size of the previous one, until an interval is found with both ends outside the slice, or until a pre-determined limit is reached. The parameters p (an integer, which determines the limit of slice size) and w (the estimate of typical slice width) are hyperparameters.

To sample from the interval \(I = (L, R)\), such that the sample \(x\) satisfies \(y < f(x)\), we use the “Shrinkage” procedure, which reduces the size of the interval after rejecting a trial point, as defined in [1] Fig. 5. pp.716. This algorithm consists of the following steps:

1. \(\bar{L} = L\) and \(\bar{R} = R\)

   2. Repeat:

      1. \(U \sim uniform(0, 1)\)
      2. \(x_1 = \bar{L} + U (\bar{R} - \bar{L})\)
      3. if \(y < f(x_1)\) and \(Accept(x_1)\), exit loop else: if \(x_1 < x_0\), then \(\bar{L} = x_1\) else \(\bar{R} = x_1\)

Intuitively, we uniformly sample a trial point from the interval I, and subsequently shrink the interval each time a trial point is rejected.

The Accept(x_1) check is required to guarantee detailed balance. We shall refer to this check as the Acceptance Check. Intuitively, it tests whether starting the doubling expansion at x_1 leads to an earlier termination compared to starting it from the current state x_0. The procedure works backward through the intervals that the doubling expansion would pass through to arrive at I when starting from x_1, checking that none of them has both ends outside the slice. The algorithm is described in [1] Fig. 6. pp.717 and it consists of the following steps:

1. \(\hat{L} = L\) and \(\hat{R} = R\) and \(D = False\)

   2. while \(\hat{R} - \hat{L} > 1.1w\):

      1. M = \((\hat{L} + \hat{R})/2\)
      2. if {\(x_0 < M\) and \(x_1 >= M\)} or {\(x_0 >= M\) and :math:` x_1 < M`}, then \(D = True\)
      3. if \(x_1 < M\), then \(\hat{R} = M\) else, \(\hat{L} = M\)
      4. if \(D\) and \(y >= f(\hat{L})\) and \(y >= f(\hat{R})\), then reject proposal

   3. If the proposal is not rejected in the previous loop, accept it

The multiplication by 1.1 in the while condition in Step 2 guards against possible round-off errors. The variable D tracks whether the intervals that would be generated from x_1 differ from those leading to x_0: when they don’t, time is saved by omitting the subsequent check.

To avoid floating-point underflow, we implement the suggestion advanced in [1] pp.712. We use the log pdf of the un-normalised posterior (\(g(x) = log(f(x))\)) instead of \(f(x)\). In doing so, we use an auxiliary variable \(z = log(y) = g(x0) - \epsilon\), where \(\epsilon \sim \text{exp}(1)\) and define the slice as \(S = {x : z < g(x)}\).

Extends SingleChainMCMC.

References

[1]	Neal, R.M., 2003. Slice sampling. The annals of statistics, 31(3), pp.705-767. https://doi.org/10.1214/aos/1056562461

ask()

See SingleChainMCMC.ask().

current_slice_height()

Returns current height value used to define the current slice.

expansion_steps()

Returns integer used for limiting interval expansion.

in_initial_phase()

For methods that need an initial phase (see needs_initial_phase()), this method returns True if the method is currently configured to be in its initial phase. For other methods a NotImplementedError is returned.

n_hyper_parameters()

See TunableMethod.n_hyper_parameters().

name()

See pints.MCMCSampler.name().

needs_initial_phase()

Returns True if this method needs an initial phase, for example an adaptation-free period for adaptive covariance methods, or a warm-up phase for DREAM.

needs_sensitivities()

Returns True if this methods needs sensitivities to be passed in to tell along with the evaluated logpdf.

replace(current, current_log_pdf, proposed=None)

Replaces the internal current position, current LogPDF, and proposed point (if any) by the user-specified values.

This method can only be used once the initial position and LogPDF have been set (so after at least 1 round of ask-and-tell).

This is an optional method, and some samplers may not support it.

set_expansion_steps(p)

Set integer for limiting interval expansion.

set_hyper_parameters(x)

The hyper-parameter vector is [width, expansion steps].

See TunableMethod.set_hyper_parameters().

set_initial_phase(in_initial_phase)

For methods that need an initial phase (see needs_initial_phase()), this method toggles the initial phase algorithm. For other methods a NotImplementedError is returned.

set_width(w)

Sets the width for generating the interval.

This can either be a single number or an array with the same number of elements as the number of variables to update.

tell(fx)

See pints.SingleChainMCMC.tell().

width()

Returns the width used for generating the interval.