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

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

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ℝ₊))

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 = LogDensityRejectErrors(ADgradient(:ForwardDiff, P));

LogDensityRejectErrors{InvalidLogDensityException,LogDensityProblems.ForwardDiffLogDensity{TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:β, :σ),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.TransformTuple{NamedTuple{(:β, :σ),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.TransformTuple{NamedTuple{(:β, :σ),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);

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([154.68144259166908, 6.1879685245233675, -0.0335400815421743, 1.5747045382544886], -1089.3699060031731, 3, DynamicHMC.DoubledTurn, 0.9956208592189933, 7), NUTS_Transition{Array{Float64,1},Float64}([155.06336625059575, 5.775236246745802, -0.34737977466365993, 1.6503933809162654], -1091.3701813334249, 3, DynamicHMC.DoubledTurn, 0.9116578504808491, 7), NUTS_Transition{Array{Float64,1},Float64}([154.98838747908349, 5.723563791147264, -0.2581881761383544, 1.6703005616874123], -1089.4414592986238, 2, DynamicHMC.AdjacentTurn, 0.9761295123241143, 7), NUTS_Transition{Array{Float64,1},Float64}([154.93329241870524, 5.7449871662394925, -0.4475666981488615, 1.7056799868587578], -1092.763196570722, 2, DynamicHMC.DoubledTurn, 0.842098026483578, 3), NUTS_Transition{Array{Float64,1},Float64}([154.6396883506302, 5.877063049463691, -0.20218484166364487, 1.696825291936204], -1092.458008487654, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([154.68698458652275, 5.820397962714695, -0.14203863661287164, 1.670289818170209], -1090.2776641427254, 3, DynamicHMC.DoubledTurn, 0.9598060436145998, 7), NUTS_Transition{Array{Float64,1},Float64}([154.54599069921792, 5.927105166725718, 0.07084436300542657, 1.6119598390156338], -1090.1222029166204, 1, DynamicHMC.AdjacentTurn, 0.8721741083525206, 3), NUTS_Transition{Array{Float64,1},Float64}([154.49080237294976, 6.024977120821876, 0.2507069066886138, 1.6552703534406723], -1090.2285415096237, 1, DynamicHMC.AdjacentTurn, 0.7930806699665314, 3), NUTS_Transition{Array{Float64,1},Float64}([154.58691853768957, 6.021386444496842, 0.2582109349202482, 1.6794485936658554], -1090.3725067097685, 2, DynamicHMC.DoubledTurn, 0.9130508431518797, 3), NUTS_Transition{Array{Float64,1},Float64}([154.80209435140958, 5.876279235940016, 0.17656167759237656, 1.600652497329514], -1090.568775831649, 2, DynamicHMC.AdjacentTurn, 0.99120793876878, 7)  …  NUTS_Transition{Array{Float64,1},Float64}([154.10451580574897, 5.930719893213091, 0.06996837892880195, 1.6585981256332705], -1093.2631000552976, 3, DynamicHMC.DoubledTurn, 0.798666202954589, 7), NUTS_Transition{Array{Float64,1},Float64}([154.19851587297907, 5.483320976521171, 0.10093012795253782, 1.6874211971123527], -1091.312750407877, 2, DynamicHMC.DoubledTurn, 0.8875220769555429, 3), NUTS_Transition{Array{Float64,1},Float64}([154.0916716222212, 5.401279683267812, 0.3587295050422174, 1.6842727823891492], -1092.7074067688022, 2, DynamicHMC.DoubledTurn, 0.91022937237642, 3), NUTS_Transition{Array{Float64,1},Float64}([155.51843012587662, 6.110425031342413, -0.46760489158069707, 1.6566851502073183], -1092.1864449251714, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([154.94430243281627, 6.186019942953062, -0.3223800435858167, 1.5902439678760907], -1092.2564857320162, 2, DynamicHMC.AdjacentTurn, 0.9902381179551946, 7), NUTS_Transition{Array{Float64,1},Float64}([154.17786873206566, 5.571160067758131, 0.3483212006691731, 1.625645830095357], -1091.6069020360947, 3, DynamicHMC.DoubledTurn, 0.9840680824130674, 7), NUTS_Transition{Array{Float64,1},Float64}([155.54115372552806, 6.487786270211303, -0.3377394907775441, 1.6679493754321295], -1094.6726163307173, 2, DynamicHMC.AdjacentTurn, 0.740801550868251, 7), NUTS_Transition{Array{Float64,1},Float64}([154.70965096034342, 6.34620215303118, -0.1760675766977283, 1.6848082635970425], -1094.8748243981952, 2, DynamicHMC.AdjacentTurn, 0.9826149831996952, 7), NUTS_Transition{Array{Float64,1},Float64}([154.60712042398876, 5.327718597513535, -0.0031178985921916647, 1.6117986067265087], -1093.9363211096552, 3, DynamicHMC.AdjacentTurn, 0.8532649285464331, 11), NUTS_Transition{Array{Float64,1},Float64}([154.45006897794553, 6.1877664061505575, 0.5333578863359139, 1.6387100224581679], -1094.407617590381, 3, DynamicHMC.DoubledTurn, 0.8047057764971309, 7)], NUTS sampler in 4 dimensions
  stepsize (ϵ) ≈ 0.648
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.357463491409066, 0.28007702228690057, 0.22060922821385556, 0.047842895509854544]
)

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.68144259166908, 6.1879685245233675, -0.0335400815421743], σ = 4.829314529595091)
 (β = [155.06336625059575, 5.775236246745802, -0.34737977466365993], σ = 5.209028556614351)
 (β = [154.98838747908349, 5.723563791147264, -0.2581881761383544], σ = 5.313764671265668) 
 (β = [154.93329241870524, 5.7449871662394925, -0.4475666981488615], σ = 5.505127809329699)
 (β = [154.6396883506302, 5.877063049463691, -0.20218484166364487], σ = 5.456596762965159) 

Extract the parameter posterior means: β,

posterior_β = mean(first, posterior)

3-element Array{Float64,1}:
 154.61847426754244    
   5.8299075628152215  
  -0.016146952236143312

then σ:

posterior_σ = mean(last, posterior)

5.095700613912254

Effective sample sizes (of untransformed draws)

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

1×4 Array{Float64,2}:
 1000.0  1000.0  682.411  1000.0

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.61847426754244, 5.8299075628152215, -0.016146952236143312]
 5.095700613912254                                               

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