$ ipython notebook --pylab="inline"
import numpy
"""
A complicated function consisting of 3 2D gaussians.
"""
c_x1 = 0.5; c_y1 = 0.2; w_x1 = 0.15; w_y1 = 0.04
c_x2 = 0.7; c_y2 = 0.7; w_x2 = 0.05; w_y2 = 0.05
c_x3 = 0.3; c_y3 = 0.7; w_x3 = 0.05; w_y3 = 0.07
def gaussian(x, y):
    return \
    exp( -( ((c_x1 - x)/w_x1)**2. + ((c_y1 - y)/w_y1)**2. )/2 ) + \
    exp( -( ((c_x2 - x)/w_x2)**2. + ((c_y2 - y)/w_y2)**2. )/2 ) + \
    exp( -( ((c_x3 - x)/w_x3)**2. + ((c_y3 - y)/w_y3)**2. )/2 )

resolution = 50.
inds = indices((resolution,resolution))/resolution
imshow( gaussian(*inds).T, origin="lower", interpolation="nearest", extent=(0,1,0,1) );

"""
A very simple Metropolis-Hastings random walker in 2D
"""
def random_walker(func, n_steps=1000, start=[0,0], burn_in=100):
    points = []
    point = start
    for i in arange(n_steps):
        # Our jumping distribution is the [-0.1,0.1] Uniform distribution around the current point
        new_point = point + (rand(2)-0.5)/5.
        accept = func(*new_point) / func(*point)
        if accept >= 1 or rand() <= accept:
            point = new_point
            if i > burn_in: points.append(point)
    return asarray(points)

"""
A function to plot the points
"""
def plot_points(points):
    h2,_,_ = histogram2d(points[:,0], points[:,1], [resolution,resolution], range=[[0,1],[0,1]])
    imshow(h2.T, origin="lower", extent=(0,1,0,1), interpolation='nearest')

# Yes, you're up!
plot_points(random_walker(gaussian, 1e4))

# Time it
%time random_walker(gaussian, 1e5);

# Fill out the missing pieces
starts = rand(10,2)
points = []
for start in starts:
    #points.append( ... )
    points.append( random_walker(gaussian, 1e4, start=start) )
plot_points( concatenate(points) )
$ ipcluster start -n 4
from IPython.parallel import Client
c = Client()

%px %pylab

%%px
start = rand(2)
#points = ...
points = random_walker(gaussian, 1e4, start=start)

plot_points(concatenate( c[:]["points"] ))
$ sudo easy_install Cython
%load_ext cythonmagic

%%px
%%cython
cimport cython
from math import pow, exp
"""
A complicated function consisting of 3 2D gaussians, now in C with Cython.
"""
cdef double c_x1 = 0.5, c_y1 = 0.2, w_x1 = 0.15, w_y1 = 0.04
cdef double c_x2 = 0.7, c_y2 = 0.7, w_x2 = 0.05, w_y2 = 0.05
cdef double c_x3 = 0.3, c_y3 = 0.7, w_x3 = 0.05, w_y3 = 0.07
def gaussian_cython(x, y):
    return \
    exp( -( pow((c_x1 - x)/w_x1,2) + pow((c_y1 - y)/w_y1,2)/2 )) + \
    exp( -( pow((c_x2 - x)/w_x2,2) + pow((c_y2 - y)/w_y2,2)/2 )) + \
    exp( -( pow((c_x3 - x)/w_x3,2) + pow((c_y3 - y)/w_y3,2)/2 ))


# Your timing code here

%%px
start = rand(2)
points = random_walker(gaussian_cython, 1e5, start=start)

plot_points(concatenate( c[:]["points"] ))

# Go for it!

import cloud

def start_walker(start):
    return random_walker(gaussian, 1e5, start=start)

starts = rand(10,2)
jids = cloud.map( start_walker, starts )

plot_points(concatenate(cloud.result(jids)))