The gibbs sampler is an iterative conditional sampler from multidimensional probability density functions (PDFs). The resulting sample is plotted as a scatter plot with the Matplotlib module.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
#Author: Flavio C. Coelho from math import * from RandomArray import * from matplotlib.pylab import * n=10000 rho=.99 #correlation #Means m1 = 10 m2 = 20 #Standard deviations s1 = 1 s2 = 1 #Initialize vectors x=zeros(n, Float) y=zeros(n, Float) sd=sqrt(1-rho**2) # the core of the method: sample recursively from two normal distributions # Tthe mean for the current sample, is updated at each step. for i in range(1,n): x[i] = normal(m1+rho*(y[i-1]-m2)/s2,s1*sd) y[i] = normal(m2+rho*(x[i-1]-m1)/s1,s2*sd) scatter(x,y,marker='d',c='r') title('Amostrador de Gibbs') xlabel('x') ylabel('y') grid() show()
This little algorithm is essencial for stochastic simulations where you need to make inferences from complex multidimensional probability density function which cannot be integrated algebraically.
Since its an iterative method it might be speeded up with Pyrex or inlining the loop in C with weave.
If anyone would like to try it please post the implementation.