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

cgibbs.pyx

'''
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+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 + 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

'''
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 + 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 +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)