Polynomial weight model model

using DynamicHMCModels

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

Import the dataset.

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

544 rows × 4 columns

	height	weight	age	male
	Float64⍰	Float64⍰	Float64⍰	Int64⍰
1	151.765	47.8256	63.0	1
2	139.7	36.4858	63.0	0
3	136.525	31.8648	65.0	0
4	156.845	53.0419	41.0	1
5	145.415	41.2769	51.0	0
6	163.83	62.9926	35.0	1
7	149.225	38.2435	32.0	0
8	168.91	55.48	27.0	1
9	147.955	34.8699	19.0	0
10	165.1	54.4877	54.0	1
11	154.305	49.8951	47.0	0
12	151.13	41.2202	66.0	1
13	144.78	36.0322	73.0	0
14	149.9	47.7	20.0	0
15	150.495	33.8493	65.3	0
16	163.195	48.5627	36.0	1
17	157.48	42.3258	44.0	1
18	143.942	38.3569	31.0	0
19	121.92	19.6179	12.0	1
20	105.41	13.948	8.0	0
21	86.36	10.4893	6.5	0
22	161.29	48.9879	39.0	1
23	156.21	42.7227	29.0	0
24	129.54	23.5868	13.0	1
⋮	⋮	⋮	⋮	⋮

Use only adults and standardize

df2 = filter(row -> row[:age] >= 18, df);
df2[:weight] = convert(Vector{Float64}, df2[:weight]);
df2[:weight_s] = (df2[:weight] .- mean(df2[:weight])) / std(df2[:weight]);
df2[:weight_s2] = df2[:weight_s] .^ 2;

352-element Array{Float64,1}:
 0.19280615097192907
 1.7349763044783346
 4.1325600023146265
 1.5549757699350777
 0.3308042660851645
 7.7736359966100865
 1.0919440933750901
 2.6392838898982456
 2.45691571770544
 2.163583954566629
 ⋮
 0.2005950450601446
 0.748125332136797
 0.37244351134337544
 0.416550377673148
 0.0999553970190397
 2.7690646673659187
 1.2340428738476732
 1.9741712709264059
 1.3641165169515888

Show the first six rows of the dataset.

first(df2, 6)

6 rows × 6 columns

	height	weight	age	male	weight_s	weight_s2
	Float64⍰	Float64	Float64⍰	Int64⍰	Float64	Float64
1	151.765	47.8256	63.0	1	0.439097	0.192806
2	139.7	36.4858	63.0	0	-1.31718	1.73498
3	136.525	31.8648	65.0	0	-2.03287	4.13256
4	156.845	53.0419	41.0	1	1.24699	1.55498
5	145.415	41.2769	51.0	0	-0.575156	0.330804
6	163.83	62.9926	35.0	1	2.78812	7.77364

Then define a structure to hold the data: observables, covariates, and the degrees of freedom for the prior.

Linear regression model $y ∼ Xβ + ϵ$, where $ϵ ∼ N(0, σ²)$ IID.

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

Then make the type callable with the parameters as a single argument.

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

Setup data and inits.

N = size(df2, 1)
X = hcat(ones(N), hcat(df2[:weight_s], df2[:weight_s2]));
y = convert(Vector{Float64}, df2[:height])
p = ConstraintHeightProblem(y, X);
p((β = [1.0, 2.0, 3.0], σ = 1.0))

-4.0013242576346947e6

Use a function to return the transformation (as it varies with the number of covariates).

problem_transformation(p::ConstraintHeightProblem) =
    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-m4.5d.ConstraintHeightProblem{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-m4.5d.ConstraintHeightProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,4,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-m4.5d.ConstraintHeightProblem{Array{Float64,1},Array{Float64,2}}}},Float64},Float64,4},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 4, w/ chunk size 4)

Draw samples.

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

MCMC, adapting ϵ (75 steps)
0.002 s/step ...done
MCMC, adapting ϵ (25 steps)
0.0014 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00055 s/step ...done
MCMC, adapting ϵ (100 steps)
0.00045 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00031 s/step ...done
MCMC, adapting ϵ (400 steps)
0.00026 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00022 s/step ...done
MCMC (1000 steps)
0.00029 s/step ...done
(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([154.91044555979818, 5.846332268039342, -0.11793249252761165, 1.5938368385631814], -1089.2033900872514, 2, DynamicHMC.DoubledTurn, 0.9779704823919669, 3), NUTS_Transition{Array{Float64,1},Float64}([155.05877892173837, 5.866283401884533, -0.10373443170623664, 1.6627807441663542], -1090.2759611051933, 2, DynamicHMC.DoubledTurn, 0.8226797160507567, 3), NUTS_Transition{Array{Float64,1},Float64}([155.170252038849, 5.847998542864838, -0.1720343085345211, 1.678620081499957], -1089.7050799930835, 1, DynamicHMC.DoubledTurn, 0.8815672606451286, 1), NUTS_Transition{Array{Float64,1},Float64}([153.98543969147758, 5.945424422179758, 0.4613947480334161, 1.606507909629199], -1092.0791344673605, 4, DynamicHMC.DoubledTurn, 0.9455597232844019, 15), NUTS_Transition{Array{Float64,1},Float64}([154.48945952043996, 5.904720103664172, -0.11218493037323546, 1.5821669497446305], -1090.8001691560523, 3, DynamicHMC.AdjacentTurn, 0.997664881951924, 15), NUTS_Transition{Array{Float64,1},Float64}([154.795120581829, 6.0349500863878855, 0.04168395246961183, 1.6590567405941943], -1089.390753438194, 3, DynamicHMC.AdjacentTurn, 0.9674406773520347, 11), NUTS_Transition{Array{Float64,1},Float64}([154.8719831660404, 5.935284604342271, -0.25629076101830045, 1.6082061757552215], -1089.8452861024723, 2, DynamicHMC.DoubledTurn, 0.9198377784787191, 3), NUTS_Transition{Array{Float64,1},Float64}([154.93200286811464, 6.159583800099805, 0.00442791110598012, 1.6329683223066234], -1089.765144452292, 3, DynamicHMC.DoubledTurn, 0.9233387352512202, 7), NUTS_Transition{Array{Float64,1},Float64}([155.0290541889373, 5.679611811574531, -0.31508426589984384, 1.6187958070850441], -1089.9219478751252, 3, DynamicHMC.DoubledTurn, 0.9926608831478537, 7), NUTS_Transition{Array{Float64,1},Float64}([154.2323836459341, 6.028517725433524, 0.24014881829926565, 1.677566446502208], -1091.549843068341, 3, DynamicHMC.AdjacentTurn, 0.8238216421706424, 15)  …  NUTS_Transition{Array{Float64,1},Float64}([154.3813908755673, 6.1148929412601145, 0.23207666731468063, 1.6179343540976492], -1090.3864449215464, 2, DynamicHMC.DoubledTurn, 0.9241672579234756, 3), NUTS_Transition{Array{Float64,1},Float64}([154.60793815271728, 6.143908072490823, 0.08551862424213413, 1.5757436346545894], -1089.9583184618355, 3, DynamicHMC.AdjacentTurn, 0.9254830211474497, 11), NUTS_Transition{Array{Float64,1},Float64}([154.35664920591634, 6.155921159424847, 0.08888142275465613, 1.67294437927654], -1090.0566016943724, 2, DynamicHMC.AdjacentTurn, 0.9783194459435247, 7), NUTS_Transition{Array{Float64,1},Float64}([154.8275851398653, 5.838665583135128, -0.1754772863595826, 1.6004887382650814], -1089.2310887299434, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([154.51098940642353, 5.882210892025473, 0.10880666330366029, 1.6569794202821877], -1088.3603902183243, 3, DynamicHMC.AdjacentTurn, 0.9625902559491621, 11), NUTS_Transition{Array{Float64,1},Float64}([155.0054904482057, 5.583784706385072, -0.1064088364031319, 1.6488173556552919], -1092.6097227649202, 2, DynamicHMC.DoubledTurn, 0.6064259214087931, 3), NUTS_Transition{Array{Float64,1},Float64}([154.85535990782603, 5.917289670921763, -0.09209302135768907, 1.6255526956210837], -1090.398651772277, 2, DynamicHMC.DoubledTurn, 0.8784962823634345, 3), NUTS_Transition{Array{Float64,1},Float64}([154.40367101357086, 5.64826644476097, 0.09915722292401635, 1.6382370457195847], -1089.5838920347705, 3, DynamicHMC.DoubledTurn, 0.8763526262360447, 7), NUTS_Transition{Array{Float64,1},Float64}([154.32564340963143, 5.971852331776687, 0.0530348440952014, 1.6133003885959294], -1089.203467654542, 3, DynamicHMC.DoubledTurn, 0.9107294433261979, 7), NUTS_Transition{Array{Float64,1},Float64}([154.8765118926338, 5.795610604465332, -0.10619144297939508, 1.64597464630313], -1088.7020545480125, 3, DynamicHMC.DoubledTurn, 0.96178022401109, 7)], NUTS sampler in 4 dimensions
  stepsize (ϵ) ≈ 0.72
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.3752025367163858, 0.2633843119411405, 0.2420178367010307, 0.04788192864997219]
)

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}:
 (β = [154.91044555979818, 5.846332268039342, -0.11793249252761165], σ = 4.922599962380172)
 (β = [155.05877892173837, 5.866283401884533, -0.10373443170623664], σ = 5.273955995168341)
 (β = [155.170252038849, 5.847998542864838, -0.1720343085345211], σ = 5.358157047304406)
 (β = [153.98543969147758, 5.945424422179758, 0.4613947480334161], σ = 4.9853714273204055)
 (β = [154.48945952043996, 5.904720103664172, -0.11218493037323546], σ = 4.865487663368754)

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

3-element Array{Float64,1}:
 154.6033832764677
   5.8369981954555685
  -0.015692240670415133

then σ:

posterior_σ = mean(last, posterior)

5.0920312345723255

Effective sample sizes (of untransformed draws)

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

NUTS-specific statistics

NUTS_statistics(chain)

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 154.609019750 0.36158389 0.0057171433 0.0071845548 1000
   b1   5.838431778 0.27920926 0.0044146860 0.0048693502 1000
   b2  -0.009985954 0.22897191 0.0036203637 0.0047224478 1000
sigma   5.110136300 0.19096315 0.0030193925 0.0030728192 1000

Quantiles:
          2.5%        25.0%        50.0%       75.0%        97.5%
    a 153.92392500 154.3567500 154.60700000 154.8502500 155.32100000
   b1   5.27846200   5.6493250   5.83991000   6.0276275   6.39728200
   b2  -0.45954687  -0.1668285  -0.01382935   0.1423620   0.43600905
sigma   4.76114350   4.9816850   5.10326000   5.2300450   5.51500975
";

"\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 154.609019750 0.36158389 0.0057171433 0.0071845548 1000\n   b1   5.838431778 0.27920926 0.0044146860 0.0048693502 1000\n   b2  -0.009985954 0.22897191 0.0036203637 0.0047224478 1000\nsigma   5.110136300 0.19096315 0.0030193925 0.0030728192 1000\n\nQuantiles:\n          2.5%        25.0%        50.0%       75.0%        97.5%\n    a 153.92392500 154.3567500 154.60700000 154.8502500 155.32100000\n   b1   5.27846200   5.6493250   5.83991000   6.0276275   6.39728200\n   b2  -0.45954687  -0.1668285  -0.01382935   0.1423620   0.43600905\nsigma   4.76114350   4.9816850   5.10326000   5.2300450   5.51500975\n"

Extract the parameter posterior means: β,

[posterior_β, posterior_σ]

2-element Array{Any,1}:
  [154.6033832764677, 5.8369981954555685, -0.015692240670415133]
 5.0920312345723255

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