#!/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()
[...] Check it out! [...]