ariddell · January 25, 2013 02:19 · May 11, 2012 · May 11, 2012 · May 11, 2012 · May 11, 2012
diff --git a/simplegibbs_cython.pyx b/simplegibbs_cython.pyx
@@ -16,10 +16,8 @@ cdef extern from "math.h":
     double sqrt(double)
 
 cdef extern from "gsl/gsl_rng.h":
-    ctypedef struct gsl_rng_type:
-        pass
-    ctypedef struct gsl_rng:
-        pass
+    ctypedef struct gsl_rng_type
+    ctypedef struct gsl_rng
 
     gsl_rng_type *gsl_rng_mt19937
     gsl_rng *gsl_rng_alloc(gsl_rng_type * T) nogil

diff --git a/simplegibbs_cython.pyx b/simplegibbs_cython.pyx
@@ -14,9 +14,6 @@ from numpy cimport *
 
 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:

diff --git a/simplegibbs_cython.pyx b/simplegibbs_cython.pyx
@@ -8,30 +8,32 @@ using conditional distributions:
 x|y \sim Gamma(3, y^2 +4)
 y|x \sim Normal(\frac{1}{1+x}, \frac{1}{2(1+x)})
 '''
-cimport cpython
+cimport cython
 import numpy as np
 from numpy cimport *
 
 cdef extern from "math.h":
     double sqrt(double)
 
-cdef double Sqrt(double n):
-    return sqrt(n)
+#cdef double Sqrt(double n):
+#    return sqrt(n)
 
 cdef extern from "gsl/gsl_rng.h":
-    ctypedef struct gsl_rng_type:   
+    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)
+    gsl_rng *gsl_rng_alloc(gsl_rng_type * T) nogil
 
-cdef gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)
-
 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)
+
+@cython.boundscheck(False)
 def gibbs(int N=20000,int thin=500):
     cdef: 
         double x=0
@@ -43,6 +45,7 @@ def gibbs(int N=20000,int thin=500):
     for i from 0 <= i < N:
         for j from 0 <= j < thin:
             x = gamma(r,3,1.0/(y*y+4))
-            y = gaussian(r,1.0/Sqrt(x+1))
-        samples[i] = x,y
-    return samples
+            y = gaussian(r,1.0/sqrt(x+1))
+        samples[i,0] = x
+        samples[i,1] = y
+    return samples
diff --git a/simplegibbs.pyx → simplegibbs_cython.pyx b/simplegibbs.pyx → simplegibbs_cython.pyx
diff --git a/simplegibbs.pyx b/simplegibbs.pyx
@@ -0,0 +1,48 @@
+'''
+Gibbs sampler for function:
+
+f(x,y) = x x^2 \exp(-xy^2 - y^2 + 2y - 4x)
+
+using conditional distributions:
+
+x|y \sim Gamma(3, y^2 +4)
+y|x \sim Normal(\frac{1}{1+x}, \frac{1}{2(1+x)})
+'''
+cimport cpython
+import numpy as np
+from numpy cimport *
+
+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 gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)
+
+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)
+
+def gibbs(int N=20000,int thin=500):
+    cdef: 
+        double x=0
+        double y=0
+        int i, j
+        ndarray[float64_t, ndim=2] samples
+
+    samples = np.empty((N,thin))
+    for i from 0 <= i < N:
+        for j from 0 <= j < thin:
+            x = gamma(r,3,1.0/(y*y+4))
+            y = gaussian(r,1.0/Sqrt(x+1))
+        samples[i] = x,y
+    return samples