Slice Sampling - Doubling MCMC¶
- class pints.SliceDoublingMCMC(x0, sigma0=None)[source]¶
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:
Calculate the pdf (\(f(x0)\)) of the current sample (\(x0\)).
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.
Find an interval (\(I = (L, R)\)) around \(x0\) that contains all, or much, of the slice.
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:
\(U \sim uniform(0, 1)\)
\(L = x_0 - wU\)
\(R = L + w\)
\(K = p\)
- while \(K > 0\) and \({y < f(L) or y < f(R)}\):
\(V \sim uniform(0, 1)\)
if \(V < 0.5\), then \(L = L - (R - L)\) else, \(R = R + (R - L)\)
\(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 parametersp
(an integer, which determines the limit of slice size) andw
(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:
- \(\bar{L} = L\) and \(\bar{R} = R\)
- Repeat:
\(U \sim uniform(0, 1)\)
\(x_1 = \bar{L} + U (\bar{R} - \bar{L})\)
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 theAcceptance Check
. Intuitively, it tests whether starting the doubling expansion atx_1
leads to an earlier termination compared to starting it from the current statex_0
. The procedure works backward through the intervals that the doubling expansion would pass through to arrive atI
when starting fromx_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:- \(\hat{L} = L\) and \(\hat{R} = R\) and \(D = False\)
- while \(\hat{R} - \hat{L} > 1.1w\):
M = \((\hat{L} + \hat{R})/2\)
if {\(x_0 < M\) and \(x_1 >= M\)} or {\(x_0 >= M\) and :math:` x_1 < M`}, then \(D = True\)
if \(x_1 < M\), then \(\hat{R} = M\) else, \(\hat{L} = M\)
if \(D\) and \(y >= f(\hat{L})\) and \(y >= f(\hat{R})\), then reject proposal
If the proposal is not rejected in the previous loop, accept it
The multiplication by
1.1
in thewhile
condition in Step 2 guards against possible round-off errors. The variableD
tracks whether the intervals that would be generated fromx_1
differ from those leading tox_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
- in_initial_phase()¶
For methods that need an initial phase (see
needs_initial_phase()
), this method returnsTrue
if the method is currently configured to be in its initial phase. For other methods aNotImplementedError
is returned.
- 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 totell
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_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 aNotImplementedError
is returned.
- set_width(w)[source]¶
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.