## Playing around with timings of Darren Wilkinson's Gibbs sampler example¶

http://darrenjw.wordpress.com/2011/07/16/gibbs-sampler-in-various-languages-revisited/

In [147]:
import random, math, time
import scipy.stats.distributions as ssd
import numpy.random as nr
import numpy

def gibbs(N=5000, thin=100):
x=0
y=0

out = []
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(2*x+2))
out += [[x, y]]
out = numpy.array(out)
return(out)

def gibbs1(N=5000, thin=100):
x=0
y=0
out = []
for i in range(N):
for j in range(thin):
x = nr.gamma(3, 1.0 / (y * y + 4))
y = nr.normal(1.0 / (x + 1), 1.0 / math.sqrt(2 * x + 2))
out += [[x, y]]
out = numpy.array(out)
return(out)

In [136]:
%load_ext cythonmagic

The cythonmagic extension is already loaded. To reload it, use:

In [154]:
%%cython -a
# non-optimize cython hack, probably some way to reach down to numpy c libraries
cimport cython
import numpy
cimport numpy
from libc.math cimport sqrt
@cython.cdivision(True)
def gibbs_cython(int N=5000, int thin=100):
cdef double x = 0, y = 0
out = []
for i in range(N):
for j in range(thin):
x = numpy.random.gamma(3, 1.0 / (y * y + 4))
y = numpy.random.normal(1.0 / (x + 1), 1.0 / sqrt(2 * x + 2))
out += [[x, y]]
out = numpy.array(out)
return(out)

Out[154]:

Generated by Cython 0.20.1 on Mon Feb 24 23:46:37 2014

 1: # non-optimize cython hack, probably some way to reach down to numpy c libraries
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 2: cimport cython
 3: import numpy
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_numpy, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 4: cimport numpy
 5: from libc.math cimport sqrt
 6: @cython.cdivision(True)
 7: def gibbs_cython(int N=5000, int thin=100):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_1gibbs_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_1gibbs_cython = {__Pyx_NAMESTR("gibbs_cython"), (PyCFunction)__pyx_pw_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_1gibbs_cython, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(0)};
static PyObject *__pyx_pw_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_1gibbs_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
int __pyx_v_N;
int __pyx_v_thin;
PyObject *__pyx_r = 0;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("gibbs_cython (wrapper)", 0);
{
static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_N,&__pyx_n_s_thin,0};
PyObject* values[2] = {0,0};
if (unlikely(__pyx_kwds)) {
Py_ssize_t kw_args;
const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
switch (pos_args) {
case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case  0: break;
default: goto __pyx_L5_argtuple_error;
}
kw_args = PyDict_Size(__pyx_kwds);
switch (pos_args) {
case  0:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_N);
if (value) { values[0] = value; kw_args--; }
}
case  1:
if (kw_args > 0) {
PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_thin);
if (value) { values[1] = value; kw_args--; }
}
}
if (unlikely(kw_args > 0)) {
if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "gibbs_cython") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
} else {
switch (PyTuple_GET_SIZE(__pyx_args)) {
case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
case  0: break;
default: goto __pyx_L5_argtuple_error;
}
}
if (values[0]) {
__pyx_v_N = __Pyx_PyInt_As_int(values[0]); if (unlikely((__pyx_v_N == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
} else {
__pyx_v_N = ((int)5000);
}
if (values[1]) {
__pyx_v_thin = __Pyx_PyInt_As_int(values[1]); if (unlikely((__pyx_v_thin == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
} else {
__pyx_v_thin = ((int)100);
}
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
__Pyx_RaiseArgtupleInvalid("gibbs_cython", 0, 0, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_L3_error:;
__Pyx_RefNannyFinishContext();
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_r = __pyx_pf_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_gibbs_cython(__pyx_self, __pyx_v_N, __pyx_v_thin);
int __pyx_lineno = 0;
const char *__pyx_filename = NULL;
int __pyx_clineno = 0;

/* function exit code */
__Pyx_RefNannyFinishContext();
return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_gibbs_cython(CYTHON_UNUSED PyObject *__pyx_self, int __pyx_v_N, int __pyx_v_thin) {
double __pyx_v_x;
double __pyx_v_y;
PyObject *__pyx_v_out = NULL;
CYTHON_UNUSED int __pyx_v_i;
CYTHON_UNUSED int __pyx_v_j;
PyObject *__pyx_r = NULL;
__Pyx_RefNannyDeclarations
__Pyx_RefNannySetupContext("gibbs_cython", 0);
/* … */
/* function exit code */
__pyx_L1_error:;
__Pyx_XDECREF(__pyx_t_1);
__Pyx_XDECREF(__pyx_t_6);
__Pyx_XDECREF(__pyx_t_7);
__Pyx_XDECREF(__pyx_t_9);
__pyx_r = NULL;
__pyx_L0:;
__Pyx_XDECREF(__pyx_v_out);
__Pyx_XGIVEREF(__pyx_r);
__Pyx_RefNannyFinishContext();
return __pyx_r;
}
/* … */
__pyx_tuple__7 = PyTuple_Pack(7, __pyx_n_s_N, __pyx_n_s_thin, __pyx_n_s_x, __pyx_n_s_y, __pyx_n_s_out, __pyx_n_s_i, __pyx_n_s_j); if (unlikely(!__pyx_tuple__7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_tuple__7);
__Pyx_GIVEREF(__pyx_tuple__7);
/* … */
__pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f_1gibbs_cython, NULL, __pyx_n_s_cython_magic_b5f9c8fec2244fa850); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
if (PyDict_SetItem(__pyx_d, __pyx_n_s_gibbs_cython, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 7; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 8:     cdef double x = 0, y = 0
  __pyx_v_x = 0.0;
__pyx_v_y = 0.0;

 9:     out = []
  __pyx_t_1 = PyList_New(0); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_v_out = __pyx_t_1;
__pyx_t_1 = 0;

 10:     for i in range(N):
  __pyx_t_2 = __pyx_v_N;
for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_2; __pyx_t_3+=1) {
__pyx_v_i = __pyx_t_3;

 11:         for j in range(thin):
    __pyx_t_4 = __pyx_v_thin;
for (__pyx_t_5 = 0; __pyx_t_5 < __pyx_t_4; __pyx_t_5+=1) {
__pyx_v_j = __pyx_t_5;

 12:             x = numpy.random.gamma(3, 1.0 / (y * y + 4))
      __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_numpy); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_random); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_gamma); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_6 = PyFloat_FromDouble((1.0 / ((__pyx_v_y * __pyx_v_y) + 4.0))); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_7 = PyTuple_New(2); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_7);
__Pyx_INCREF(__pyx_int_3);
PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_int_3);
__Pyx_GIVEREF(__pyx_int_3);
PyTuple_SET_ITEM(__pyx_t_7, 1, __pyx_t_6);
__Pyx_GIVEREF(__pyx_t_6);
__pyx_t_6 = 0;
__pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_7, NULL); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
__pyx_t_8 = __pyx_PyFloat_AsDouble(__pyx_t_6); if (unlikely((__pyx_t_8 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_v_x = __pyx_t_8;

 13:             y = numpy.random.normal(1.0 / (x + 1), 1.0 / sqrt(2 * x + 2))
      __pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_numpy); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_7 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_random); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_7);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_7, __pyx_n_s_normal); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
__pyx_t_7 = PyFloat_FromDouble((1.0 / (__pyx_v_x + 1.0))); if (unlikely(!__pyx_t_7)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_7);
__pyx_t_1 = PyFloat_FromDouble((1.0 / sqrt(((2.0 * __pyx_v_x) + 2.0)))); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_9 = PyTuple_New(2); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_9);
PyTuple_SET_ITEM(__pyx_t_9, 0, __pyx_t_7);
__Pyx_GIVEREF(__pyx_t_7);
PyTuple_SET_ITEM(__pyx_t_9, 1, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_1);
__pyx_t_7 = 0;
__pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_9, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
__pyx_t_8 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_8 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 13; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_y = __pyx_t_8;
}

 14:         out += [[x, y]]
    __pyx_t_1 = PyFloat_FromDouble(__pyx_v_x); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_9 = PyFloat_FromDouble(__pyx_v_y); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_9);
__pyx_t_6 = PyList_New(2); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
PyList_SET_ITEM(__pyx_t_6, 0, __pyx_t_1);
__Pyx_GIVEREF(__pyx_t_1);
PyList_SET_ITEM(__pyx_t_6, 1, __pyx_t_9);
__Pyx_GIVEREF(__pyx_t_9);
__pyx_t_1 = 0;
__pyx_t_9 = 0;
__pyx_t_9 = PyList_New(1); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_9);
PyList_SET_ITEM(__pyx_t_9, 0, __pyx_t_6);
__Pyx_GIVEREF(__pyx_t_6);
__pyx_t_6 = 0;
__pyx_t_6 = PyNumber_InPlaceAdd(__pyx_v_out, __pyx_t_9); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
__Pyx_DECREF_SET(__pyx_v_out, __pyx_t_6);
__pyx_t_6 = 0;
}

 15:     out = numpy.array(out)
  __pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_numpy); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_9 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_array); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_9);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_6 = PyTuple_New(1); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_6);
__Pyx_INCREF(__pyx_v_out);
PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_v_out);
__Pyx_GIVEREF(__pyx_v_out);
__pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_9, __pyx_t_6, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 15; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_DECREF_SET(__pyx_v_out, __pyx_t_1);
__pyx_t_1 = 0;

 16:     return(out)
  __Pyx_XDECREF(__pyx_r);
__Pyx_INCREF(__pyx_v_out);
__pyx_r = __pyx_v_out;
goto __pyx_L0;

In [155]:
def gibbs1(N=5000, thin=100):
x=0
y=0
out = []
for i in range(N):
for j in range(thin):
x = nr.gamma(3, 1.0 / (y * y + 4))
y = nr.normal(1.0 / (x + 1), 1.0 / math.sqrt(2 * x + 2))
out += [[x, y]]
out = numpy.array(out)
return(out)

In [156]:
t = time.time()
out = gibbs()
print(time.time() - t)

%timeit gibbs1()

t = time.time()
out1 = gibbs1()
print(time.time() - t)

%timeit gibbs_cython()

t = time.time()
out2 = gibbs_cython()
print(time.time() - t)

2.900670051574707
1 loops, best of 3: 659 ms per loop
0.6385161876678467
1 loops, best of 3: 262 ms per loop
0.29333996772766113

In [157]:
from pylab import *
figure(1)
ax = subplot(111)
df = numpy.hstack([out, out1, out2])
df = pandas.DataFrame(df)
p = df.plot(ax=ax, kind='density', alpha=0.5, linewidth=3, style=['-', '--', 'x'] * 2)

In [ ]: