Linear regression Gibbs sampler

Specific Domains Probabilistic programming
1

I tried to copy this Python sampler Gibbs sampling for Bayesian linear regression in Python | Kieran R Campbell - blog but my Julia implementation returns unexpected numbers.

What am I doing wrong?


using Distributions, LinearAlgebra, Random, DataFrames, UnPack



function draw_β₀(y, x, β₁, τ, μ₀, τ₀)
    precision = τ₀ + τ*n
    location = (τ₀ * μ₀ + τ * sum(y .- β₁ * x)) / precision
    rand(Normal(location, 1/sqrt(precision)))
end

function draw_β₁(y, x, β₀, τ, μ₁, τ₁)
    precision = τ₁ + τ * (x ⋅ x)
    location = (τ₁ * μ₁ + τ * sum((y .- β₀) .* x)) / precision
    rand(Normal(location, 1/sqrt(precision)))
end

function draw_τ(y, x, β₀, β₁, α, β)
    α = α + n / 2
    resid = y .- β₀ .- β₁ * x
    β = β + (resid ⋅ resid) / 2
    rand(Gamma(α, β)) # alpha beta parametrization
end


function gibbs(y, x, iters, θ₀, hypers)
    @unpack β₀, β₁, τ =  θ₀
    @unpack μ₀, τ₀, μ₁, τ₁, α, β = hypers
    trace = DataFrame(; β₀, β₁, τ)
    for _ in 1:iters
        β₀ = draw_β₀(y, x, β₁, τ, μ₀, τ₀)
        β₁ = draw_β₁(y, x, β₀, τ, μ₁, τ₁)
        τ = draw_τ(y, x, β₀, β₁, α, β)
        push!(trace, (; β₀, β₁, τ))
    end
    return trace
end


Random.seed!(1)
n = 50
Xq = [ones(n) 4rand(n)]
βq = [-10., 5.]
τq = 2.
ϵq = sqrt(τq) * randn(n)
μq = Xq*βq
zq = μq + ϵq
iters = 500
θ₀ = (β₀=0., β₁=0., τ=2.)
hypers = (μ₀=0., τ₀=1., μ₁=0., τ₁=1., α=2., β=1.)


trace = gibbs(zq, Xq, iters, θ₀, hypers)


mean.(eachcol(trace))
3-element Array{Float64,1}:
 -8.183514118517033e71
  4.0211361378904415e71
  2.34582711781041e148
2

Hi @jzr ,
both, in numpy.random, as well as Distributions.jl, the Gamma distribution is parametrized by shape and scale parameters (not rate), as far as I can tell.
In your function draw_τ, β seems to play the role of a rate parameter (as in the original implementation).
However the original implementation calls np.random.gamma with 1/β, you call Gamma with β.

(I didn’t run your code, so not sure whether changing that fixes it)

2 Likes
3

Yeah, compare with Seven Lines of Julia (examples sought) - #26 by mschauer

2 Likes