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

%%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)

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)

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)