Slice Sampler (2D)

#!/usr/bin/python
"""
Author: Jeremy M. Stober
Program: SLICE.PY
Date: Sunday, August 17 2008
Description: Slice sample method as described in "Monte Carlo Statistical Methods" by Robert and Casella.
"""
 
import os, sys, getopt, pdb, string
from numpy import *
import pylab
 
def slice2d(niters, func, finv, x):
    samples = [x]
    for i in range(niters):
        u = random.uniform(low=0.0, high=func(x))    
        x = random.uniform(low=0.0, high=finv(u))
        samples.append(x)
    return samples
 
def func(x):
    return 0.5 * exp(-sqrt(x))
 
def finv(u):
    return log(2 * u)**2
 
 
comparison = [func(x) for x in arange(0,30,0.1)]
samples = slice2d(5000, func,finv, 0)
pylab.cla()
pylab.hist(samples,bins=100,normed=True)
pylab.plot(arange(0,30,0.1), comparison, '-', color='red')
pylab.xlim(0,30)
pylab.show()
Share and Enjoy:
  • Digg
  • del.icio.us
  • Facebook
  • Google Bookmarks
  • Reddit
  • Technorati
  • Furl
  • StumbleUpon
  • Tumblr
  • TwitThis
One Response to “Slice Sampler (2D)”
  1.  
Leave a Reply