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