Research Memorandum No. 1209

01/06/2020

Unbounded Slice Sampling
Mochihashi Daichi (The Institute of Statistical Mathematics)

The Institute of Statistical Mathematics
10-3, Midori-cho, Tachikawa,Tokyo, 190-8562, Japan

Unbounded Slice Sampling
Daichi Mochihashi
The Institute of Statistical Mathematics
daichi@ism.ac.jp
2020-1-5 (Sun)
Abstract
Slice sampling is an eﬃcient Markov Chain Monte Carlo algorithm to sample from anunnormalized density with acceptance ratio always 1. However, when the variable to sample
is unbounded, its “stepping-out” heuristic works only locally, making it diﬃcult to uniformly
explore possible candidates. This paper proposes a simple change-of-variable method to slice
sample an unbounded variable equivalently from [ 0, 1).

1Introduction
Slice sampling [1] is one of the Markov Chain Monte Carlo (MCMC) methods to sample

from one-dimensional distribution. Because its acceptance probability is always 1 as opposed tofamous Metropolis-Hastings algorithm [2], slice sampling is widely employed in modern statistics
and machine learning problems, especially for sampling hyperparameters of a statistical model
where each hyperparameter is often univariate but has a nontrivial likelihood surface.
When the variable to sample is not bounded, an algorithm called “stepping-out” [3] is usuallyemployed to adaptively adjust the interval to sample from based on the current value of thevariable. However, this “stepping-out” works only locally, thus will be potentially trapped to
local modes in multimodal likelihoods. In contrast, in this paper we propose a simple method touniformly sample from an unbounded variable through an appropriate reparametrization with aunit interval [ 0, 1).

2

Slice sampling
Slice sampling iteratively draws x from a probability distributionf (x)(1)
Z!
without knowing the normalizing constant Z = f (x)dx by the following general algorithm,p(x) =starting from an initial value x(0) :Algorithm 1 General slice sampling
1: for t = 1 · · · T do
2:ℓ = f (x(t−1) )3:ρ ∼ Uniform[0, ℓ)

4:x(t) ∼ Uniform({x : f (x) > ρ})

5: end for1

(a) Sigmoid map p =1
1 + e−x(b) Logit map x =p1−pFigure 1: One-to-one map between p ∈ [ 0, 1) and unbounded x.
In practice, Step 4 of the Algorithm 1 is not trivial, and usually “stepping-out” heuristic [1]is employed to determine the interval from which x is uniformly drawn, starting from an initialinterval enclosing x that is adaptively expanded.
However, the eﬃciency of “stepping-out” depends on setting the initial interval whose appropriate scale for the problem is not known in advance. Moreover, it works only locally: if f (x)has multiple and distant modes, it tends to wander around a single mode and will spend a very
long time to explore diﬀerent modes that might have higher probabilities.
On the other hand, if x is bounded and known to always reside in the interval [st, ed ], slice
sampling is very easy using the following algorithm:
Algorithm 2 Bounded slice sampling
1: for t = 1 · · · T doℓ = f (x(t−1) )3:ρ ∼ Uniform[0, ℓ)
4:while true do
5:x ∼ Uniform[st, ed ]6:
if f (x) > ρ then7:x(t) = x; break
8:
else
9:
if x < x(t−1) then10:
st := x
11:
else
12:
ed := x13:
end if14:
end if15:
end while
16: end for2:/* generate a candidate *//* modify interval */This is an eﬃcient binary search that quickly draws x uniformly over the region [st, ed ] wheref (x) > ρ.
2

3

Unbounded slice sampling

Algorithm 2 is eﬀective for a bounded x that is known to reside in [st, ed ]. However, evenwhen x is unbounded, we can sample x through an one-to-one map between [ 0, 1) and R.
General unbounded variate For example, using a sigmoid map
1
1 + e−x"#1x = − log−1pp=i.e.

(2)(3)

for p ∈ [ 0, 1) shown in Figure 1(a), for each x ∈ R we can associate p ∈ [ 0, 1) and vice versa.
Since

dx
1=,dp
p(1 − p)we can alternatively sample p in place of x using the exchange of variables:
""
##dx
1
1f (x)dx = f (x) dp = f − log−1dpdppp(1 − p)

(4)

(5)

Using (5) for the Algorithm 2 in place of f (x) can sample p ∈ [ 0, 1), from which we can recoverx by the relationship (3).
Positive variate

When x > 0, for example we can exploit a mapx=thereforep1−p(6)dx1=,dp(1 − p)2

leadingdx
f (x)dx = f (x) dp = fdp"p1−pto sample p ∈ [ 0, 1) and obtain x by the relationship (6).

(7)
#1dp
(1 − p)2(8)

Using this reparameterization, we can slice sample an unbounded x essentially by a “proba-bilistic binary search” always on [ 0, 1). We included sample of MATLAB and C codes for ourunbounded slice sampling in Appendix B and C.
In practice, sometimes we cannot compute a corresponding p for very large (or small) x using
equation (2). Therefore, instead of (2) we can use
p = σ(x/A) =1
1 + e−x/A(9)for some constant A > 0 as shown in Figure 2. While every x ∈ R is mathematically equivalent,

notice that usually the parameter to sample from a statistical model roughly1 resides in thisregime; even if x will exceed over this regime, we can safely rescale x to fit into a decent interval.
We empirically confirmed that A = 100 usually works fine for −1000 < x < 1000 (see Section 4).

4

Experiments
1 This is important, because we do not enforce a strict interval to sample x in this paper.

3

500

1500

400
1000

300
200
500

100
0-101

2

30
250

4

300

350

400

450

500

550

2(a) f (x) = exp(−x(x−1)(x−2)(x−3.5))(b) f (x) = exp(−(x−500) /10)

10000
8000
6000
4000
2000
00

200

400

600

800

1000

1200

2(c) f (x) = exp(−(x−1000) /100)
Figure 3: Unbounded slice sampling of x ∈ R.
We conducted some experiments to confirm that
our unbounded slice sampler correctly samples froma given distribution. Using the C code shown in Appendix A, we can easily sample from p(x).
For a general unbounded case x ∈ R, Figure 3 shows

the plots of sampled x from (a) f (x) = exp(−(x −500)2 /10), (b) f (x) = exp(−(x − 500)2 /10), and (c)f (x) = exp(−(x − 1000)2 /100). Note that p(x) will

be high on diﬀerent regimes for each case: roughly
0 < x < 4 for (a), x ≈ 500 for (b), and x ≈ 1000 for(c), while these regimes are not known in advance,
and they have diﬀerent variances. Figure 3 clearlyshows we can correctly sample from p(x), even whenx is very large using the map (9) with A = 100. The

Figure 2: Modified sigmoid map p =1.
−x/1001+eaverage step of “binary search”, i.e. the number of function evaluations for the loop in Step 4 of
the Algorithm 2, is 11.44 for (a), 16.48 for (b), and 9.34 for (c).
For a positive case x ∈ R+ , Figure 4 shows the results for (a) f (x) = x4 e−x and (b) f (x) =exp(−(x−1000)2 /100); (a) means p(x) ∼ Gamma(5, 1).
Comparison with stepping-out

To confirm the advantage of our unbounded slice sampling,

we also conducted an experiment to sample from a multimodal distribution. Specifically, Fig-4350

6000

300

5000

250
4000

200
3000

150
2000

100
1000

50000

5104 −x(a) f (x) = x e15200200

400

600

800

1000

12002(= Gamma(5, 1))(b) f (x) = exp(−(x−1000) /100)

Figure 4: Unbounded slice sampling of x > 0.
1500

8006001000400500
2000-505100-515(a) Unbounded slice sampling

051015(b) Stepping-out algorithm (width=1)

Figure 5: Comparison with “stepping-out” on a Gaussian mixture model p(x) = 0.8N (x|0, 12 ) +0.2N (x|10, 12 ) using 10, 000 samples.
ure 5 shows a draw of 10, 000 samples from a Gaussian mixture model p(x) = 0.8N (x|0, 12 ) +0.2N (x|10, 12 ) using unbounded slice sampler and “stepping-out” algorithm, respectively.
Clearly, “stepping-out” is stuck into the first mode because the sampling starts from x = 1,and this tendency remained the same over multiple simulations. Moreover, because “steppingout” linearly increases the interval to sample from, sampling from a distant density shown in
Figure 3(c) required about 2, 000 function evaluations at first to reach the mode, while unboundedslice sampling required only 9.34 as described before, thanks to the eﬃcient “probabilistic binarysearch” on R.

References[1] Radford M. Neal. Slice sampling. Annals of Statistics, pages 705–741, 2003.
[2] W. R. Gilks, S. Richardson, and D. J. Spiegelhalter. Markov Chain Monte Carlo in Practice.
Chapman & Hall / CRC, 1996.
[3] David J. C. MacKay. Information Theory, Inference, and Learning Algorithms. Cambridge
University Press, 2003.

5

Appendix
A

C code of easy slice sampling#include "slice.h"doublef (double x, void *arg){
return - x * (x - 1) * (x - 2) * (x - 3.5);}
main () {int i, N = atoi (argv[1]);double x = 1;
for (i = 0; i < N; i++) {x = slice_sample (x, f, NULL);
printf("%lf\n", x);}}B

C code for unbounded slice sampling/**slice.c*Unbounded slice sampling with ease.
*$Id: unbounded.tex,v 1.17 2020/01/03 01:23:08 daichi Exp $**/#include <stdio.h>#include <stdlib.h>#include <math.h>#include "random.h"static double A = 100.0;
typedef double (*likfun)(double x, void *args);
static double
expand (double p)/* p -> x */{return - A * log (1 / p - 1);}
static double
shrink (double x)/* x -> p */{return 1 / (1 + exp(- x / A));}static double
expandp (double p)/* p -> (x > 0) */{return p / (1 - p);}static double
shrinkp (double x)/* (x > 0) -> p */{return x / (1 + x);

6}
doubleslice_sample (double x, likfun loglik, void *arg){
double st = 0, ed = 1;
double r, rnew, slice, newlik;
int iter, maxiter = 1000;r = shrink (x);slice = (*loglik)(x, arg) - log (A * r * (1 - r)) + log (RANDOM);
for (iter = 0; iter < maxiter; iter++){rnew = unif (st, ed);newlik = (*loglik)(expand(rnew), arg)
- log (A * rnew * (1 - rnew));
if (newlik > slice)
return expand (rnew);else if (rnew > r)
ed = rnew;
else if (rnew < r)
st = rnew;
else
return x;}
fprintf(stderr, "slice_sample: max iteration %d reached.\n", maxiter);return x;}
doubleslice_sample_positive (double x, likfun loglik, void *arg){
double st = 0, ed = 1;
double r, rnew, slice, newlik;
int iter, maxiter = 1000;r = shrinkp (x);slice = (*loglik)(x, arg) - 2 * log (1 - r) + log (RANDOM);
for (iter = 0; iter < maxiter; iter++){rnew = unif (st, ed);newlik = (*loglik)(expandp(rnew), arg) - 2 * log (1 - rnew);if (newlik > slice)
return expandp (rnew);else if (rnew > r)
ed = rnew;
else if (rnew < r)
st = rnew;
else
return x;}
fprintf(stderr, "slice_sample: max iteration %d reached.\n", maxiter);return x;}#if 1
doublef (double x, void *arg){
return - x * (x - 1) * (x - 2) * (x - 3.5);}
double7g (double x, void *arg){
return - 10 * (x + 1000) * (x + 1000);}
doubleh (double p, void *arg){
double a = 2, b = 3;
return (a - 1) * log (p) + (b - 1) * log (1 - p);}intmain (int argc, char *argv[]){
int i, N = atoi(argv[1]);double x = 0.5;
for (i = 0; i < N; i++){x = slice_sample (x, g, NULL);// x = slice_sample_positive (x, f, NULL);
printf("%.4lf\n", x);}}#endif

CMATLAB code for unbounded slice sampling

function xnew = slice_sample (x,likfun,varargin)% xnew = slice_sample (x,likfun,varargin)% Unbounded slice sampling in MATLAB.
% $Id: unbounded.tex,v 1.17 2020/01/03 01:23:08 daichi Exp $global A;
A = 100;
st = 0; ed = 1;
maxiter = 1000;% function bodyr = shrink (x);slice = feval (likfun,x,varargin{:}) - log (A * r * (r - 1)) + log (rand());% slice samplingfor iter = 1:maxiter
rnew = unif (st,ed);xnew = expand (rnew);newlik = feval (likfun,xnew,varargin{:}) - log (A * rnew * (rnew - 1));
if (newlik > slice)
return;elseif (rnew > r)
ed = rnew;
elseif (rnew < r)
st = rnew;
else
return;end
end
fprintf(stderr,’slice_sample: max iteration %d reached\n’,maxiter);return;function p = shrink(x) % x -> pglobal A;
p = 1 / (1 + exp (- x / A));
function x = expand(p) % p -> xglobal A;
x = - A * log (1 / p - 1);

8

function xnew = slice_sample_positive (x,likfun,varargin)% xnew = slice_sample_positive (x,likfun,varargin)% Unbounded slice sampling in MATLAB. (only positive)% $Id: unbounded.tex,v 1.17 2020/01/03 01:23:08 daichi Exp $st = 0; ed = 1;
maxiter = 1000;% function bodyr = shrinkp (x);slice = feval (likfun,x,varargin{:}) - 2 * log (1 - r) + log (rand());% slice samplingfor iter = 1:maxiter
rnew = unif (st,ed);xnew = expandp (rnew);newlik = feval (likfun,xnew,varargin{:}) - 2 * log (1 - rnew);if (newlik > slice)
return;elseif (rnew > r)
ed = rnew;
elseif (rnew < r)
st = rnew;
else
return;end
end
fprintf(stderr,’slice_sample: max iteration %d reached\n’,maxiter);return;function p = shrinkp (x) % (x > 0) -> p
p = x / (1 + x);function x = expandp (p) % p -> (x > 0)x = p / (1 - p);9