public
Last active

Comparison of MCMC implementations in Python and Cython. This is discussed here: http://pyinsci.blogspot.com/2010/12/efficcient-mcmc-in-python.html

  • Download Gist
cgibbs.pyx
Cython
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
'''
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()
cgibbs1.pyx
Cython
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
'''
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()
cgibbs2.pyx
Cython
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
'''
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()
gibbs.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
'''
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)
gibbs1.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
'''
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)
gibbs2.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
'''
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)

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.