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
```