Linear regression

using DynamicHMCModels

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

Import the dataset.

snippet 5.4

wd = CSV.read(rel_path("..", "data", "WaffleDivorce.csv"), delim=';')
df = convert(DataFrame, wd);

mean_ma = mean(df[!, :Marriage])
df[!, :Marriage_s] = convert(Vector{Float64},
  (df[!, :Marriage]) .- mean_ma)/std(df[!, :Marriage]);

mean_mam = mean(df[!, :MedianAgeMarriage])
df[!, :MedianAgeMarriage_s] = convert(Vector{Float64},
  (df[!, :MedianAgeMarriage]) .- mean_mam)/std(df[!, :MedianAgeMarriage]);

50-element Array{Float64,1}:
 -0.6062895051354262 
 -0.6866992538271283 
 -0.20424076167692148
 -1.4103869920524357 
  0.5998567252400879 
 -0.28465051036862354
  1.2431347147736962 
  0.43903722785668664
  2.9317394372994143 
  0.27821773047328247
  ⋮                  
 -0.6866992538271283 
 -0.6866992538271283 
 -2.214484478969445  
  0.6802664739317872 
  0.27821773047328247
 -0.12383101298522224
 -0.8475187512105296 
  0.19780798178158324
 -1.4907967407441378 

Show the first six rows of the dataset.

first(df[!, [1, 7, 14,15]], 6)

6 rows × 4 columns

Model $y ∼ Xβ + ϵ$, where $ϵ ∼ N(0, σ²)$ IID. Student prior on σ

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

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

function (problem::m_5_3)(θ)
    @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(Normal(0, 1), X[3]) # b1 = X[3]
    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[!, :Marriage_s], df[!, :MedianAgeMarriage_s]);
y = convert(Vector{Float64}, df[!, :Divorce])
p = m_5_3(y, X);
p((β = [1.0, 2.0, 3.0], σ = 1.0))

-2222.175273500088

Write a function to return properly dimensioned transformation.

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

problem_transformation (generic function with 1 method)

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

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

ForwardDiff AD wrapper for TransformedLogDensity of dimension 4, w/ chunk size 4

Tune and sample.

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

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([9.44859839780487, -0.38822720992007287, -1.2063909833808442, 0.3060722742054358], -101.5016375107753, 2, DynamicHMC.AdjacentTurn, 0.9555681429989432, 7), NUTS_Transition{Array{Float64,1},Float64}([9.477748270509286, -0.22179098622197196, -0.9739309388564404, 0.33493555952395154], -100.18687471355072, 3, DynamicHMC.DoubledTurn, 0.9676507339397684, 7), NUTS_Transition{Array{Float64,1},Float64}([9.660700097092928, -0.25905977771917593, -1.3664411707177266, 0.4699387900407548], -99.34451345923559, 2, DynamicHMC.DoubledTurn, 0.9858957869169793, 3), NUTS_Transition{Array{Float64,1},Float64}([9.901473245396385, -0.685280558297139, -1.4309784425536518, 0.5020299483898931], -102.54894318505188, 1, DynamicHMC.AdjacentTurn, 0.7153640173271075, 3), NUTS_Transition{Array{Float64,1},Float64}([9.611805549501947, 0.10058266887514294, -1.3762127315396824, 0.39397522770119164], -101.86775403249462, 2, DynamicHMC.AdjacentTurn, 0.9194731679246563, 7), NUTS_Transition{Array{Float64,1},Float64}([9.632009348797693, -0.3597538304948612, -1.4440707920872293, 0.4504951157595679], -100.98014720078153, 2, DynamicHMC.DoubledTurn, 0.9400394114929629, 3), NUTS_Transition{Array{Float64,1},Float64}([9.497904993523461, -0.18382762383042134, -1.3608648272096973, 0.45746120028442605], -98.94934097766524, 2, DynamicHMC.AdjacentTurn, 0.9723065665259837, 7), NUTS_Transition{Array{Float64,1},Float64}([9.670926433644762, -0.45583801001007873, -1.5494999616307734, 0.30467701550787524], -99.17101173474788, 3, DynamicHMC.DoubledTurn, 0.9941309802823081, 7), NUTS_Transition{Array{Float64,1},Float64}([9.563095152771764, -0.13878242450858486, -1.0637828459963596, 0.3072849999686285], -100.41750307928406, 3, DynamicHMC.AdjacentTurn, 0.856852782875545, 15), NUTS_Transition{Array{Float64,1},Float64}([9.651755118014494, -0.04547826039031169, -0.888785739815467, 0.38380733621986746], -99.04324008323319, 2, DynamicHMC.AdjacentTurn, 0.960638503352037, 7)  …  NUTS_Transition{Array{Float64,1},Float64}([9.72871302863776, -0.5036983606124653, -1.8000975674115414, 0.37581294074594385], -99.96374627796848, 2, DynamicHMC.DoubledTurn, 0.8672844572135731, 3), NUTS_Transition{Array{Float64,1},Float64}([9.918965677864591, 0.031625860465330086, -0.6479073777168369, 0.4043617636718365], -101.05692212257783, 3, DynamicHMC.DoubledTurn, 0.9516132425613971, 7), NUTS_Transition{Array{Float64,1},Float64}([9.6211505124098, -0.4107979399844264, -1.5032186924926172, 0.22701626919251125], -104.44322872876324, 3, DynamicHMC.DoubledTurn, 0.7251810126352941, 7), NUTS_Transition{Array{Float64,1},Float64}([9.652591702022093, -0.2953741042051475, -1.1910488192167263, 0.3435665847421631], -100.07703922944597, 2, DynamicHMC.DoubledTurn, 0.9052906834599247, 3), NUTS_Transition{Array{Float64,1},Float64}([9.786421450345038, -0.1750790955522674, -1.3037080220365704, 0.41489081978325476], -99.02040830711982, 3, DynamicHMC.DoubledTurn, 0.9268681226228843, 7), NUTS_Transition{Array{Float64,1},Float64}([9.638181409247348, -0.14547320571315575, -1.4723258729910864, 0.3475958490675849], -100.4741619593947, 2, DynamicHMC.DoubledTurn, 0.815649776492292, 3), NUTS_Transition{Array{Float64,1},Float64}([9.8656736339535, -0.6513476579805358, -1.1556781639097735, 0.4921719023871601], -101.34601581611095, 2, DynamicHMC.AdjacentTurn, 0.7762830561593663, 7), NUTS_Transition{Array{Float64,1},Float64}([9.569871220959001, 0.35133403015349357, -1.144362105838841, 0.347622862857632], -101.96024528828711, 3, DynamicHMC.DoubledTurn, 0.9800892748238047, 7), NUTS_Transition{Array{Float64,1},Float64}([9.64317820093146, -0.6805983828949609, -1.586702739210891, 0.2333747775486135], -103.08230011011821, 3, DynamicHMC.DoubledTurn, 0.9865568480406027, 7), NUTS_Transition{Array{Float64,1},Float64}([9.565589171827192, 0.1969432457210854, -0.9640338591631108, 0.5028715467905426], -100.3336189876059, 3, DynamicHMC.DoubledTurn, 1.0, 7)], NUTS sampler in 4 dimensions
  stepsize (ϵ) ≈ 0.75
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.19281852893239879, 0.3063603664098533, 0.3075931116578144, 0.1021851290592024]
)

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}:
 (β = [9.44859839780487, -0.38822720992007287, -1.2063909833808442], σ = 1.3580804571869312) 
 (β = [9.477748270509286, -0.22179098622197196, -0.9739309388564404], σ = 1.3978503041810197)
 (β = [9.660700097092928, -0.25905977771917593, -1.3664411707177266], σ = 1.599896260635256) 
 (β = [9.901473245396385, -0.685280558297139, -1.4309784425536518], σ = 1.65207148902368)    
 (β = [9.611805549501947, 0.10058266887514294, -1.3762127315396824], σ = 1.4828638142742578) 

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

3-element Array{Float64,1}:
  9.675756289333941 
 -0.2170463173406174
 -1.2351324361764504

then σ:

posterior_σ = mean(last, posterior)

1.5031993651540325

Effective sample sizes (of untransformed draws)

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

1×4 Array{Float64,2}:
 864.581  904.233  875.131  886.911

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.91, min/25%/median/75%/max: 0.44 0.86 0.95 0.99 1.0
  termination: AdjacentTurn => 16% DoubledTurn => 84%
  depth: 1 => 2% 2 => 47% 3 => 50% 4 => 1% 5 => 0%

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.69137275 0.21507432 0.0034006235 0.0038501180 1000
   bA -1.12184710 0.29039965 0.0045916216 0.0053055477 1000
   bM -0.12106472 0.28705400 0.0045387223 0.0051444688 1000
sigma  1.52326545 0.16272599 0.0025729239 0.0034436330 1000

Quantiles:
         2.5%       25.0%      50.0%      75.0%       97.5%
    a  9.2694878  9.5497650  9.6906850  9.83227750 10.11643500
   bA -1.6852295 -1.3167700 -1.1254650 -0.92889225 -0.53389157
   bM -0.6889247 -0.3151695 -0.1231065  0.07218513  0.45527243
sigma  1.2421182  1.4125950  1.5107700  1.61579000  1.89891925
";

"\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.69137275 0.21507432 0.0034006235 0.0038501180 1000\n   bA -1.12184710 0.29039965 0.0045916216 0.0053055477 1000\n   bM -0.12106472 0.28705400 0.0045387223 0.0051444688 1000\nsigma  1.52326545 0.16272599 0.0025729239 0.0034436330 1000\n\nQuantiles:\n         2.5%       25.0%      50.0%      75.0%       97.5%\n    a  9.2694878  9.5497650  9.6906850  9.83227750 10.11643500\n   bA -1.6852295 -1.3167700 -1.1254650 -0.92889225 -0.53389157\n   bM -0.6889247 -0.3151695 -0.1231065  0.07218513  0.45527243\nsigma  1.2421182  1.4125950  1.5107700  1.61579000  1.89891925\n"

Extract the parameter posterior means: [β, σ],

[posterior_β, posterior_σ]

2-element Array{Any,1}:
  [9.675756289333941, -0.2170463173406174, -1.2351324361764504]
 1.5031993651540325                                            

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