Skip to content

Instantly share code, notes, and snippets.

@fonnesbeck
Created May 11, 2012 15:27
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 3 You must be signed in to fork a gist
  • Save fonnesbeck/2660449 to your computer and use it in GitHub Desktop.
Save fonnesbeck/2660449 to your computer and use it in GitHub Desktop.
Simple Gibbs sampler in cython, adapted from http://bit.ly/J3vP63
'''
Gibbs sampler for function:
f(x,y) = x x^2 \exp(-xy^2 - y^2 + 2y - 4x)
using conditional distributions:
x|y \sim Gamma(3, y^2 +4)
y|x \sim Normal(\frac{1}{1+x}, \frac{1}{2(1+x)})
'''
cimport cython
import numpy as np
from numpy cimport *
cdef extern from "math.h":
double sqrt(double)
cdef extern from "gsl/gsl_rng.h":
ctypedef struct gsl_rng_type
ctypedef struct gsl_rng
gsl_rng_type *gsl_rng_mt19937
gsl_rng *gsl_rng_alloc(gsl_rng_type * T) nogil
cdef extern from "gsl/gsl_randist.h":
double gamma "gsl_ran_gamma"(gsl_rng * r,double,double)
double gaussian "gsl_ran_gaussian"(gsl_rng * r,double)
cdef gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)
@cython.boundscheck(False)
def gibbs(int N=20000,int thin=500):
cdef:
double x=0
double y=0
int i, j
ndarray[float64_t, ndim=2] samples
samples = np.empty((N,thin))
for i from 0 <= i < N:
for j from 0 <= j < thin:
x = gamma(r,3,1.0/(y*y+4))
y = gaussian(r,1.0/sqrt(x+1))
samples[i,0] = x
samples[i,1] = y
return samples
@fonnesbeck
Copy link
Author

Timings relative to a pure Python implementation:

In [1]: import simplegibbs_cython

In [2]: timeit simplegibbs_cython.gibbs(20000, 200)
1 loops, best of 3: 513 ms per loop

In [3]: import simplegibbs

In [4]: timeit simplegibbs.gibbs(20000, 200)
1 loops, best of 3: 26.3 s per loop

Thanks to Flavio Coelho for the original cython implementation, which I modified. The performance here compares favorably to versions coded in Julia and Rcpp:

julia> sum([@elapsed JGibbs1(20000, 200) for i=1:10])/10
0.39336202144622813

These benchmarks were run on a 2011 MacBook Air under OS X 10.7.3.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment