Skip to content

Instantly share code, notes, and snippets.

@fccoelho
Created December 28, 2010 09:37
Show Gist options
  • Save fccoelho/757090 to your computer and use it in GitHub Desktop.
Save fccoelho/757090 to your computer and use it in GitHub Desktop.
Comparison of MCMC implementations in Python and Cython. This is discussed here: http://pyinsci.blogspot.com/2010/12/efficcient-mcmc-in-python.html
'''
Pure cython version
compile with:
$ cython cgibbs.pyx
$ gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.6 -o cgibbs.so cgibbs.c
then import from python shell and call main()
'''
import random,math, time
def gibbs(int N=20000,int thin=500):
cdef double x=0
cdef double y=0
cdef int i, j
samples = []
for i in range(N):
for j in range(thin):
x=random.gammavariate(3,1.0/(y*y+4))
y=random.gauss(1.0/(x+1),1.0/math.sqrt(x+1))
samples.append((x,y))
return samples
def main():
t0 = time.time()
gibbs()
print "time %s seconds"%(time.time()-t0)
#smp = gibbs()
'''
cython+numpy version
compile with:
$ cython cgibbs1.pyx
$ gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.6 -o cgibbs1.so cgibbs1.c
then import from python shell and call main()
'''
import numpy as np
cimport numpy as np
import time
def gibbs(int N=20000,int thin=500):
cdef double x=0
cdef double y=0
cdef int i, j
samples = []
for i in range(N):
for j in range(thin):
x=np.random.gamma(3,1.0/(y*y+4))
y=np.random.normal(1.0/(x+1),1.0/np.sqrt(x+1))
samples.append((x,y))
return samples
def main():
t0 = time.time()
gibbs()
print "time %s seconds"%(time.time()-t0)
#smp = gibbs()
'''
cython + GSL version
compile with:
$ cython cgibbs1.pyx
$ gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.6 -lgsl -lblas -o cgibbs2.so cgibbs2.c
then import from python shell and call main()
'''
#declaring external GSL functions to be used
cdef extern from "math.h":
double sqrt(double)
cdef double Sqrt(double n):
return sqrt(n)
cdef extern from "gsl/gsl_rng.h":
ctypedef struct gsl_rng_type:
pass
ctypedef struct gsl_rng:
pass
gsl_rng_type *gsl_rng_mt19937
gsl_rng *gsl_rng_alloc(gsl_rng_type * T)
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)
# original Cython code
def gibbs(int N=20000,int thin=500):
cdef double x=0
cdef double y=0
cdef int i, j
samples = []
#print "Iter x y"
for i in range(N):
for j in range(thin):
x = gamma(r,3,1.0/(y*y+4))
y = 1.0/(x+1)+gaussian(r,1.0/Sqrt(x+1))
samples.append([x,y])
return samples
def main():
import time
t0 = time.time()
smp = gibbs()
print "time: %s seconds"%(time.time()-t0)
#smp = gibbs()
'''
Pure Python Version
'''
import random,math,time
def gibbs(N=20000,thin=500):
x=0
y=0
samples = []
for i in range(N):
for j in range(thin):
x=random.gammavariate(3,1.0/(y*y+4))
y=random.gauss(1.0/(x+1),1.0/math.sqrt(x+1))
samples.append((x,y))
return samples
t0 = time.time()
smp = gibbs()
print "time %s seconds"%(time.time()-t0)
'''
Python + numpy
'''
import time
import numpy as np
def gibbs(N=20000,thin=500):
x=0
y=0
samples = []
for i in range(N):
for j in range(thin):
x=np.random.gamma(3,1.0/(y*y+4))
y=np.random.normal(1.0/(x+1),1.0/np.sqrt(x+1))
samples.append((x,y))
return samples
t0 = time.time()
smp = gibbs()
print "time %s seconds"%(time.time()-t0)
'''
Python +numpy without using list to store results
'''
import random,math,time
import numpy as np
def gibbs(N=20000,thin=500):
x=0
y=0
samples = np.empty((20000,2), dtype=float)
for i in range(N):
for j in range(thin):
x=np.random.gamma(3,1.0/(y*y+4))
y=np.random.normal(1.0/(x+1),1.0/np.sqrt(x+1))
samples[i]=(x,y)
return samples
t0 = time.time()
smp = gibbs()
print "time %s seconds"%(time.time()-t0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment