5  Linear Regression

We can further apply Bayesian statistics to help with estimating the parameters of a linear regression model. Here, I’m going to use Julia’s Turing.jl package to fit the model.

using Distributions
using Plots
using Turing
using LinearAlgebra
using Random
using StatsPlots
using GLM

Random.seed!(0408)

[ Info: Precompiling GLM [38e38edf-8417-5370-95a0-9cbb8c7f171a]

TaskLocalRNG()

5.1 Generate Fake Data

I’ll generate some fake data to use as an example.

n = 1000

𝐗 = randn(n, 3)

β = [1.0, 2.0, 3.0]

f(x) = 0.5 .+ x * β

ϵ = rand(Normal(0, 1.0), n)

y = f(𝐗) + ϵ;

Here’s what we’re modeling:

\(Y_i \sim N(\alpha + \beta^T \bf X_i, \sigma^2)\)

In words, \(Y_i\) is normally distributed with an expected value of the output of a typical linear model and a variance of \(\sigma^2\).

5.2 Define Model

Next, we can use Turing’s @model macro to define our model. One important part of defining our model is setting priors.

In this model, I’m going to use pretty weak priors, and I’m going to use these same priors for all of the variables.

@model function lin_reg(x, y)
    n_feat = size(x, 2)

    #priors
    α ~ Normal(0, 2)
    b ~ MvNormal(zeros(n_feat), 3.0 * I)
    σ ~ Exponential(1)

    #estimate the likelihood
    for i ∈ eachindex(y)
        y[i] ~ Normal(α + x[i, :]' * b, σ)
    end

    #alternatively:
    # μ = α .+ x * b
    # return y ~ MvNormal(μ, σ * I)
end

lin_reg (generic function with 2 methods)

In math notation, our model is:

\(Y_i | \alpha, \beta, \sigma \sim N(\alpha + \bf X_i ' \beta, \sigma^2)\)

\(\alpha \sim N(0, 2)\)

\(\beta \sim N(\mu, \Sigma)\)

where \(\mu\) is a length-3 zero-vector and \(\Sigma\) is a diagonal matrix with 3s on the diagonal

\(\sigma \sim Exp(1)\)

5.3 Simulate the Posterior

Ok so now we have this model set up, which is what a fairly typical linear regression model might look like, and now we need to simulate the posterior. We can do that in Turing as follows:

model = lin_reg(𝐗, y)

chn = sample(model, NUTS(), MCMCThreads(), 5_000, 2);

┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC C:\Users\eric_ekholm\.julia\packages\AbstractMCMC\Es490\src\sample.jl:307
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.003125

And we can check the trace plot and the posterior distributions:

plot(chn)