Source code for qikify.controllers.slicesample

import numpy as np

[docs]def slicesample(x0, nsamples, pdf, width = 10, maxiter = 200): """Loosely based on slicesample() from MATLAB. """ dim = np.size(x0) rnd = np.zeros((nsamples,dim)) e = np.random.exponential(1,nsamples) # needed for the vertical position of the slice. RW = np.random.rand(nsamples,dim) # factors of randomizing the width RD = np.random.rand(nsamples,dim) # uniformly draw the point within the slice for i in xrange(nsamples): # A vertical level is drawn uniformly from (0,f(x0)) and used to define # the horizontal "slice". z = logpdf(x0, pdf) - e[i] # An interval [xl, xr] of width w is randomly position around x0 and then # expanded in steps of size w until both size are outside the slice. r = width * RW[i,:] xl = x0 - r xr = xl + width iteration = 0 # step out procedure is performed only when univariate samples are drawn. if dim==1: # step out to the left. while inside(xl,z, pdf) and iteration < maxiter: xl -= width iteration += 1 # step out to the right iteration = 0 while inside(xr,z, pdf) and iteration < maxiter: xr += width iteration += 1 # A new point is found by picking uniformly from the interval [xl, xr]. xp = RD[i,:]*(xr-xl) + xl # shrink the interval (or hyper-rectangle) if a point outside the # density is drawn. iteration = 0 while outside(xp,z, pdf) and iteration < maxiter: rshrink = (xp > x0) lshrink = ~rshrink xr[rshrink] = xp[rshrink] xl[lshrink] = xp[lshrink] xp = (np.random.rand(1,dim) * (xr-xl))[0] + xl # draw again iteration += 1 rnd[i,:] = x0 = xp # update the current value return (rnd[0,:] if nsamples == 1 else rnd)
[docs]def logpdf(x, pdf): fx = pdf(x) return np.log(fx) if fx > 0 else -np.inf
[docs]def inside(x,th, pdf): return logpdf(x, pdf) > th
[docs]def outside(x,th, pdf): return logpdf(x, pdf) <= th