Slice sampling algorithm

A single transition $(x,u) \rightarrow (x',u')$ of a one-dimensional slice sampling algorithm has the following steps.

(1). evaluate $P^* (x)$

(2). draw a vertical coordinate $u' \sim$ Uniform$(0,P^* (x))$

(3). create a horizontal interval $(x_l, x_r)$ enclosing $x$

    3a. draw $r \sim$ Uniform$(0,1)$

    3b. $x_l := x-rw$

    3c. $x_r := x+(1-r)w$

    3d. while $(P^* (x_l) > u')$ ${x_l := x-rw}$

    3e. while $(P^* (x_r) > u')$ ${x_r:= x+w}$

(4). loop {

(5).     draw $x' \sim$ Uniform$(x_l,x_r)$

(6).     evaluate $P^*(x')$

(7).     if $P^*(x') >u'$ break out of loop 4-9

(8).     else modify the interval $(x_l,x_r)$

    8a.     if $(x' > x){x_r:=x'}$

    8b.     else ${x_l:=x'}$

(9). }

I find the following exercise is quite simple and very good to understand how the slice sampling algorithm works. This is the exercise 29.10 of the book Information Theory, Inference, and Learning Algorithms of David Mackay.

Exercise 29.10

Question

Investigate the properties of slice sampling applied to the density function $P^* (x)$.

$$ P^* (x)= \{ \begin{array}{cc} 10 ~~~~~~as ~0\leq x<1\\\
~~0 ~~~~~~as ~1\leq x \leq 11\
\end{array} $$

$x$ is a real variable between 0.0 and 11.0. How long does it take typically for slice sampling to get from an $x$ in the peak region $x \in (0,~ 1)$ to an $x$ in the tail region $x \in (1,~11)$, and vice versa? Confirm that the probabilities of these transitions do yield an asymptotic probability density that is correct.

Solution

import numpy as np
import matplotlib.pyplot as plt
import random
import timeit

Let us see how the target probability distribution looks like.

# probability density

def Prob(x):
    if x < 0:
        return 0
    elif x < 1:
        return 10
    elif x <11 :
        return 1
    else :
        return 0

# probability density plot

x = np.arange(-2.0, 13.0, 0.1)
P=[Prob(t) for t in x]
y = np.array(P)
plt.plot(x, y)
plt.show()

png

Give the initial sample between $x_0 = 0$ and $1$, and decide the size of the step. It should be not too big nor too small. I chose it 0.5.


x0 = random.random()

# decide the size of the step. Not too small, not too big.
w = 0.5

The following is the main part of the slice sampling. Compare with the above algorithm.

r = random.random()
u = Prob(x0)*random.random()

xl = x0 - r * w

while Prob(xl) > u:
    xl = xl - w


xr = x0 + (1-r)*w

while Prob(xr) > u:
    xr = xr + w

Put them together for the step of the process.

def process(x0, w):


    # decide the size of the step. Not too small, not too big.
    w = 0.5

    u = Prob(x0)*random.random()

    cond = True
    r = random.random()
    xl = x0-r*w
    xr = x0 + (1-r)*w

    while cond:
        while Prob(xl) > u:
            xl = xl - w

        while Prob(xr) > u:
            xr = xr + w

        x1= random.uniform(xl,xr)

        if Prob(x1) > u :
            cond = False
        else :
            if x1>x0:
                xr = x1
            else:
                xl=x1
    return x1

Use the process method, I want to find the number of turns the sample comes out from the range, (0,1) to out of the range.

def numofturns(x0, w):
    turn = 0
    while x0 < 1:
        turn = turn + 1
        x0 = process(x0, 0.5)

    return turn


list =[]
for i in range(10000):
    list = list+ [numofturns(random.random(),0.5)]


plt.plot(list)
plt.show()

png

np.average(list)

10.9325

It takes around 11 truns.

Then, how about opposite? From out of the range to inside of the range, (0,1).

def numofturns2(x0, w):
    turn = 0
    while x0 > 1:
        turn = turn + 1
        x0 = process(x0, 0.5)

    return turn


list2 =[]
for i in range(10000):
    list2 = list2 + [numofturns2(random.uniform(1,11),0.5)]

plt.plot(list2)
plt.show()

png

np.average(list2)

11.1015

It takes about 11 turns as well. The difference is about 1.5 percentages. When you calculate the areas of the plot below (0,1 ), and out of the range. They match. It supports our analysis that it takes the same time for the samples to go in or out of the range, (0,1).