## A Julia version of Darren Wilkinson's R code
function JGibbs1(N::Int, thin::Int)
    mat = Array(Float64, (N, 2))
    x   = 0.
    y   = 0.
    for i = 1:N
        for j = 1:thin
            x = randg(3)/(y*y + 4)
            y = 1/(x + 1) + randn(1)[1]/sqrt(2(x + 1))
        end
        mat[i,:] = [x,y]
    end
    mat
end

load("extras/Rmath.jl")

function JGibbs2(N::Int, thin::Int)
    mat = Array(Float64, (N, 2))
    x   = 0.
    y   = 0.
    for i = 1:N
        for j = 1:thin
            x = rgamma(1,3,(y*y + 4))[1]
            y = rnorm(1, 1/(x+1),1/sqrt(2(x + 1)))[1]
        end
        mat[i,:] = [x,y]
    end
    mat
end

dJGibbs1(N::Int, thin::Int) = convert(Array{Float64,2}, darray((T,d,da)->JGibbs1(d[1],thin), Float64, (N, 2), 1))

dJGibbs2(N::Int, thin::Int) = darray((T,d,da)->JGibbs1(d[1],thin), Float64, (N, 2), 1)