This is an example from Dr. Doug Bates's slides Julia for R Programmers.
The task is to create a Gibbs sampler for the density \[ f(x, y) = k x^2 exp(-xy^2 - y^2 + 2y - 4x), x > 0 \] using the conditional distributions \[ \begin{eqnarray*} X | Y &\sim& \Gamma \left(3, \frac{1}{y^2 + 4}\right) \\ Y | X &\sim& N \left(\frac{1}{1+x}, \frac{1}{2(1+x)} \right). \end{eqnarray*} \]
This is a Julia function for the simple Gibbs sampler
using Distributions
function jgibbs(N::Integer, thin::Integer)
mat = Array(Float64, (N, 2))
x = y = 0.0
for i = 1:N
for j in 1:thin
x = rand(Gamma(3.0, 1.0 / (y * y + 4.0)))
y = rand(Normal(1.0 / (x + 1.0), 1.0 / sqrt(2.0(x + 1.0))))
end
mat[i, 1] = x
mat[i, 2] = y
end
mat
end
Generate a bivariate sample of size 10,000 with a thinning of 500. How long does it take?
#jgibbs(100, 5); # warm-up
@elapsed jgibbs(10000, 500)
Distributed version in Julia.
djgibbs(N::Integer, thin::Integer) = DArray(I -> jgibbs(map(length, I)[1], thin), (N, 2))
This function creates distributed array (different parts are stored and accessed on different processors) using an anonymous function, which is the first argument and uses the -> operator to create a function.
@elapsed djgibbs(10000, 500)
I didn't figure out how to run multi-threading in IJulia and somehow we don't see gain from the distributed code. When run in terminal mode, I see about 3 fold speedup in djgibbs.
Show system information.
versioninfo()