slice.py

def slice_sample(likelihood, x0):
    old_lik = likelihood(x0)
    old_alpha = x0
    lnt = old_lik - np.random.exponential(1)
    w = old_alpha/32.
    L = max(0, old_alpha - w*random.random())
    R = L + w
    K = 4
    while K > 0 and (lnt < likelihood(L) or lnt < likelihood(R)):
        V = random.random()
        if V < 0.5:
            if L-(R-L) < 0:
                print "L would be", L-(R-L), "R is", R
            L = max(0, L-(R-L))
        else:
            R = R+(R-L)
        K = K-1
    rej = True
    while rej:
        U = random.random()
        x1 = L+U*(R-L)
        rr = likelihood(x1)
        if lnt < rr:
            break
        else:
            if x1 < old_alpha:
                L = x1
            else:
                R = x1
    return x1