# Linear regression Gibbs sampler

**URL:** <https://discourse.julialang.org/t/linear-regression-gibbs-sampler/62817>\
**Category:** Probabilistic Programming\
**Created:** [June 12, 2021, 10:57pm UTC](https://discourse.julialang.org/t/linear-regression-gibbs-sampler/62817 "2021-06-12T22:57:31Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jzr](https://avatars.discourse-cdn.com/v4/letter/j/eb9ed0/32.png) [@jzr](https://discourse.julialang.org/u/jzr)\
**Post date:** [June 12, 2021, 10:57pm UTC](https://discourse.julialang.org/t/linear-regression-gibbs-sampler/62817/1 "2021-06-12T22:57:31Z")

</div>

I tried to copy this Python sampler [Gibbs sampling for Bayesian linear regression in Python | Kieran R Campbell - blog](https://kieranrcampbell.github.io/blog/2016/05/15/gibbs-sampling-bayesian-linear-regression.html) but my Julia implementation returns unexpected numbers.

What am I doing wrong?

```julia

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

```

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [July 3, 2021, 10:40am UTC](https://discourse.julialang.org/t/linear-regression-gibbs-sampler/62817/2 "2021-07-03T10:40:24Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![mschauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mschauer/32/13946_2.png) [@mschauer](https://discourse.julialang.org/u/mschauer)\
**Post date:** [July 3, 2021, 11:06am UTC](https://discourse.julialang.org/t/linear-regression-gibbs-sampler/62817/3 "2021-07-03T11:06:53Z")

</div>

Yeah, compare with [Seven Lines of Julia (examples sought) - #26 by mschauer](https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/26)