# Linear regression Gibbs sampler

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

``````

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

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

2 Likes