Linear regression

using DynamicHMCModels, MCMCChains

ProjDir = rel_path_d("..", "scripts", "05")
cd(ProjDir)

Import the dataset.

snippet 5.1

wd = CSV.read(rel_path("..", "data", "WaffleDivorce.csv"), delim=';')
df = convert(DataFrame, wd);
mean_ma = mean(df[:MedianAgeMarriage])
df[:MedianAgeMarriage_s] = convert(Vector{Float64},
  (df[:MedianAgeMarriage]) .- mean_ma)/std(df[:MedianAgeMarriage]);

Show the first six rows of the dataset.

first(df, 6)

6 rows × 14 columns

	Location	Loc	Population	MedianAgeMarriage	Marriage	Marriage SE	Divorce	Divorce SE	WaffleHouses	South	Slaves1860	Population1860	PropSlaves1860	MedianAgeMarriage_s
	String⍰	String⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Float64⍰	Float64
1	Alabama	AL	4.78	25.3	20.2	1.27	12.7	0.79	128	1	435080	964201	0.45	-0.60629
2	Alaska	AK	0.71	25.2	26.0	2.93	12.5	2.05	0	0	0	0	0.0	-0.686699
3	Arizona	AZ	6.33	25.8	20.3	0.98	10.8	0.74	18	0	0	0	0.0	-0.204241
4	Arkansas	AR	2.92	24.3	26.4	1.7	13.5	1.22	41	1	111115	435450	0.26	-1.41039
5	California	CA	37.25	26.8	19.1	0.39	8.0	0.24	0	0	0	379994	0.0	0.599857
6	Colorado	CO	5.03	25.7	23.5	1.24	11.6	0.94	11	0	0	34277	0.0	-0.284651

Model $y ∼ Normal(y - Xβ, σ)$. Flat prior for β, half-T for σ.

struct WaffleDivorceProblem{TY <: AbstractVector, TX <: AbstractMatrix}
    "Observations."
    y::TY
    "Covariates"
    X::TX
end

Make the type callable with the parameters as a single argument.

function (problem::WaffleDivorceProblem)(θ)
    @unpack y, X, = problem   # extract the data
    @unpack β, σ = θ            # works on the named tuple too
    ll = 0.0
    ll += logpdf(Normal(10, 10), X[1]) # a = X[1]
    ll += logpdf(Normal(0, 1), X[2]) # b1 = X[2]
    ll += logpdf(TDist(1.0), σ)
    ll += loglikelihood(Normal(0, σ), y .- X*β)
    ll
end

Instantiate the model with data and inits.

N = size(df, 1)
X = hcat(ones(N), df[:MedianAgeMarriage_s]);
y = convert(Vector{Float64}, df[:Divorce])
p = WaffleDivorceProblem(y, X);
p((β = [1.0, 2.0], σ = 1.0))

-2225.6614871340917

Write a function to return properly dimensioned transformation.

problem_transformation(p::WaffleDivorceProblem) =
    as((β = as(Array, size(p.X, 2)), σ = asℝ₊))

Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(problem_transformation(p), p)
∇P = LogDensityRejectErrors(ADgradient(:ForwardDiff, P));

LogDensityRejectErrors{InvalidLogDensityException,LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.1d1.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.1d1.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.ShiftedExp{true,Float64}}},Main.ex-m5.1d1.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 3, w/ chunk size 3)

Create an array to hold 1000 samples of 3 parameters in 4 chains

a3d = create_a3d(1000, 3, 4);
trans = as( (β = as(Array, 2), σ = asℝ));

TransformVariables.TransformNamedTuple{(:β, :σ),Tuple{TransformVariables.ArrayTransform{TransformVariables.Identity,1},TransformVariables.Identity}}((TransformVariables.ArrayTransform{TransformVariables.Identity,1}(TransformVariables.Identity(), (2,)), TransformVariables.Identity()), 3)

Sample from 4 chains and store the draws in the a3d array

for j in 1:4
  chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 3000);
  posterior = TransformVariables.transform.(Ref(problem_transformation(p)),
    get_position.(chain));
  insert_chain!(a3d, j, posterior, trans);
end;

cnames = ["a", "bA", "sigma"]
chns = Chains(a3d, cnames)

MCMC, adapting ϵ (75 steps)
5.7e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
5.8e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
5.7e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
3.9e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00013 s/step ...done
MCMC, adapting ϵ (400 steps)
3.5e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
6.8e-5 s/step ...done
MCMC (3000 steps)
7.9e-5 s/step ...done
MCMC, adapting ϵ (75 steps)
0.0003 s/step ...done
MCMC, adapting ϵ (25 steps)
8.8e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
5.0e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
4.0e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
3.2e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
6.3e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
3.6e-5 s/step ...done
MCMC (3000 steps)
5.2e-5 s/step ...done
MCMC, adapting ϵ (75 steps)
4.4e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
4.9e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
5.2e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
4.0e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
8.6e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
3.4e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
4.0e-5 s/step ...done
MCMC (3000 steps)
7.6e-5 s/step ...done
MCMC, adapting ϵ (75 steps)
4.9e-5 s/step ...done
MCMC, adapting ϵ (25 steps)
4.9e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
4.7e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
0.00014 s/step ...done
MCMC, adapting ϵ (200 steps)
3.3e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
3.2e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00025 s/step ...done
MCMC (3000 steps)
6.7e-5 s/step ...done
Object of type Chains, with data of type 1000×3×4 Array{Float64,3}

Log evidence      = 0.0
Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 1000
parameters        = a, bA, sigma

parameters
        Mean    SD   Naive SE  MCSE   ESS
    a  9.6839 0.2128   0.0034 0.0035 1000
   bA -1.0859 0.2202   0.0035 0.0033 1000
sigma  1.4934 0.1634   0.0026 0.0034 1000

cmdstan result

cmdstan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
         Mean        SD       Naive SE       MCSE      ESS
    a  9.6882466 0.22179190 0.0035068378 0.0031243061 1000
   bA -1.0361742 0.21650514 0.0034232469 0.0034433245 1000
sigma  1.5180337 0.15992781 0.0025286807 0.0026279593 1000

Quantiles:
         2.5%      25.0%     50.0%      75.0%       97.5%
    a  9.253141  9.5393175  9.689585  9.84221500 10.11121000
   bA -1.454571 -1.1821025 -1.033065 -0.89366925 -0.61711705
sigma  1.241496  1.4079225  1.504790  1.61630750  1.86642750
";

"\nIterations = 1:1000\nThinning interval = 1\nChains = 1,2,3,4\nSamples per chain = 1000\n\nEmpirical Posterior Estimates:\n         Mean        SD       Naive SE       MCSE      ESS\n    a  9.6882466 0.22179190 0.0035068378 0.0031243061 1000\n   bA -1.0361742 0.21650514 0.0034232469 0.0034433245 1000\nsigma  1.5180337 0.15992781 0.0025286807 0.0026279593 1000\n\nQuantiles:\n         2.5%      25.0%     50.0%      75.0%       97.5%\n    a  9.253141  9.5393175  9.689585  9.84221500 10.11121000\n   bA -1.454571 -1.1821025 -1.033065 -0.89366925 -0.61711705\nsigma  1.241496  1.4079225  1.504790  1.61630750  1.86642750\n"

Extract the parameter posterior means: β,

describe(chns)

Log evidence      = 0.0
Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3, 4
Samples per chain = 1000
parameters        = a, bA, sigma

Empirical Posterior Estimates
─────────────────────────────────────────────
parameters
        Mean    SD   Naive SE  MCSE   ESS
    a  9.6839 0.2128   0.0034 0.0035 1000
   bA -1.0859 0.2202   0.0035 0.0033 1000
sigma  1.4934 0.1634   0.0026 0.0034 1000

Quantiles
─────────────────────────────────────────────
parameters
        2.5%   25.0%   50.0%   75.0%   97.5%
    a  8.9484  9.5475  9.6843  9.8237 10.6310
   bA -1.9209 -1.2337 -1.0845 -0.9385 -0.3098
sigma  1.0548  1.3744  1.4811  1.5921  2.2590

Plot the chains

plot(chns)

end of m4.5d.jl#- This page was generated using Literate.jl.