Linear regression

using DynamicHMCModels

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.1d.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.1d.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.1d.WaffleDivorceProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,3},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 3, w/ chunk size 3)

Tune and sample.

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000);

MCMC, adapting ϵ (75 steps)
0.00029 s/step ...done
MCMC, adapting ϵ (25 steps)
8.0e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
6.1e-5 s/step ...done
MCMC, adapting ϵ (100 steps)
0.0002 s/step ...done
MCMC, adapting ϵ (200 steps)
3.5e-5 s/step ...done
MCMC, adapting ϵ (400 steps)
4.1e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00026 s/step ...done
MCMC (1000 steps)
4.4e-5 s/step ...done
(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([10.09010240085756, -1.528698603557108, 0.5152497762281153], -100.63804600808245, 2, DynamicHMC.DoubledTurn, 0.9905736145587047, 3), NUTS_Transition{Array{Float64,1},Float64}([9.964195901484779, -1.636424657655045, 0.46450964386731075], -101.75470898705498, 2, DynamicHMC.AdjacentTurn, 0.9741803546818707, 7), NUTS_Transition{Array{Float64,1},Float64}([9.449202181263471, -0.6703298066814263, 0.30604717948385346], -102.14241593643402, 3, DynamicHMC.DoubledTurn, 0.9686393738880585, 7), NUTS_Transition{Array{Float64,1},Float64}([9.603831507023415, -0.9629484001847674, 0.4791152105059254], -99.55815900807947, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([9.560335809568635, -1.4544082369775468, 0.40925769557658015], -98.46942427384283, 3, DynamicHMC.DoubledTurn, 0.9568729950210019, 7), NUTS_Transition{Array{Float64,1},Float64}([9.345352259494554, -1.2109428275692882, 0.39751531731206735], -99.02620813241708, 2, DynamicHMC.AdjacentTurn, 0.9824555671240853, 7), NUTS_Transition{Array{Float64,1},Float64}([9.362087899319143, -0.9603236687755798, 0.4259690777283178], -98.90752992796526, 2, DynamicHMC.DoubledTurn, 0.9768308559416777, 3), NUTS_Transition{Array{Float64,1},Float64}([10.016157168130224, -1.1984316289370407, 0.3521034759820775], -98.39314149463996, 3, DynamicHMC.DoubledTurn, 0.99778736994761, 7), NUTS_Transition{Array{Float64,1},Float64}([9.750009373194624, -0.9327193649793957, 0.30481207202807276], -98.3149080938546, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([9.700608813859244, -1.2693789366713217, 0.4517276115868591], -99.19654932328196, 3, DynamicHMC.DoubledTurn, 0.8353477588129545, 7)  …  NUTS_Transition{Array{Float64,1},Float64}([9.572924478680779, -0.9713897668751263, 0.45705382133666017], -100.70407597873188, 2, DynamicHMC.DoubledTurn, 0.8080371810067429, 3), NUTS_Transition{Array{Float64,1},Float64}([9.917147086045029, -0.8607195563153708, 0.4418289281907412], -99.98906610439903, 2, DynamicHMC.DoubledTurn, 0.8308547062968761, 3), NUTS_Transition{Array{Float64,1},Float64}([9.47865471757887, -1.2696978553099214, 0.32700277168780817], -98.14710313238314, 3, DynamicHMC.DoubledTurn, 0.995554271658425, 7), NUTS_Transition{Array{Float64,1},Float64}([9.893563308557075, -0.9599621882579061, 0.41085508350252764], -100.07895283078324, 3, DynamicHMC.DoubledTurn, 0.8752144907830336, 7), NUTS_Transition{Array{Float64,1},Float64}([9.889423912519332, -0.8651527690114478, 0.4245022505350109], -97.59337784545615, 2, DynamicHMC.DoubledTurn, 0.9893264117961259, 3), NUTS_Transition{Array{Float64,1},Float64}([9.686286522186743, -1.453116842912702, 0.3012627141307152], -99.23576823535879, 3, DynamicHMC.DoubledTurn, 0.9467675796019291, 7), NUTS_Transition{Array{Float64,1},Float64}([9.542780188140501, -1.5252950707148873, 0.2975664704073571], -99.78700125142375, 2, DynamicHMC.DoubledTurn, 0.9540761320976426, 3), NUTS_Transition{Array{Float64,1},Float64}([9.357217492929625, -1.3884951041376503, 0.3095873151135703], -100.61606507127836, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([9.608421639591732, -1.3756729421095024, 0.5098957902413701], -100.12250778201513, 2, DynamicHMC.AdjacentTurn, 0.9945363240601063, 7), NUTS_Transition{Array{Float64,1},Float64}([9.730894246806406, -1.051531022352155, 0.451453567006008], -100.22546917884051, 2, DynamicHMC.AdjacentTurn, 0.923080291501212, 5)], NUTS sampler in 3 dimensions
  stepsize (ϵ) ≈ 0.653
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.22247113616576766, 0.2090183697763557, 0.10582722275414289]
)

We use the transformation to obtain the posterior from the chain.

posterior = TransformVariables.transform.(Ref(problem_transformation(p)), get_position.(chain));
posterior[1:5]

5-element Array{NamedTuple{(:β, :σ),Tuple{Array{Float64,1},Float64}},1}:
 (β = [10.09010240085756, -1.528698603557108], σ = 1.6740565891217007)
 (β = [9.964195901484779, -1.636424657655045], σ = 1.59123372633419)
 (β = [9.449202181263471, -0.6703298066814263], σ = 1.3580463769635902)
 (β = [9.603831507023415, -0.9629484001847674], σ = 1.6146451489275682)
 (β = [9.560335809568635, -1.4544082369775468], σ = 1.505699682643222)

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

2-element Array{Float64,1}:
  9.684472745270568
 -1.0894510552632788

then σ:

posterior_σ = mean(last, posterior)

1.495399588799254

Effective sample sizes (of untransformed draws)

ess = mapslices(effective_sample_size,
                get_position_matrix(chain); dims = 1)

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.94, min/25%/median/75%/max: 0.56 0.91 0.97 0.99 1.0
  termination: AdjacentTurn => 18% DoubledTurn => 82%
  depth: 1 => 4% 2 => 48% 3 => 48%

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: β,

[posterior_β[1], posterior_β[2], posterior_σ]

3-element Array{Float64,1}:
  9.684472745270568
 -1.0894510552632788
  1.495399588799254

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