http://darrenjw.wordpress.com/2011/07/16/gibbs-sampler-in-various-languages-revisited/
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)
%load_ext cythonmagic
The cythonmagic extension is already loaded. To reload it, use: %reload_ext cythonmagic
%%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)
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_AddTraceback("_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f.gibbs_cython", __pyx_clineno, __pyx_lineno, __pyx_filename); __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_AddTraceback("_cython_magic_b5f9c8fec2244fa850f9ed1f9e85434f.gibbs_cython", __pyx_clineno, __pyx_lineno, __pyx_filename); __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;
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)
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
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)