Seven Lines of Julia (examples sought)

Agree with Tamas, but it is also nice to have a couple of very short examples. Here a Gibbs sampler sampling the posterior parameters of normal observations with mean distributed a priori as N(μ0, σ0) and precision 1/σ^2 a priori as Gamma(α, β):

using Distributions, Statistics
function mc(x, iterations, μ0 = 0., τ0 = 1e-7, α = 0.0001, β = 0.0001)
    sumx, n = sum(x), length(x)
    μ, τ = sumx/n, 1/var(x)
    samples = [(μ, τ)]
    for i in 1:iterations
        μ = rand(Normal((τ0*μ0 + τ*sumx)/(τ0 + n*τ), 1/sqrt(τ0 + n*τ)))
        τ = rand(Gamma(α + n/2, 1/(β + 0.5*sum((xᵢ-μ)^2 for xᵢ in x))))
        push!(samples, (μ, τ))
    end
    samples
end

Usage

samples = mc(rand(Normal(0.2, 1.7^(-0.5)), 1_000), 10_000)
μs, τs = first.(samples), last.(samples)
println("mean μ = ", mean(μs), " ± ", std(μs))
println("precision τ = ", mean(τs), " ± ", std(τs))

Output:

mean μ = 0.2076960556449643 ± 0.023693641548604788
precision τ = 1.8301418451485263 ± 0.08225999887566306

9 Likes