Heights problem

We estimate simple linear regression model with a half-T prior.

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);

352 rows × 4 columns

Show the first six rows of the dataset.

first(df2, 6)

6 rows × 4 columns

Half-T for σ, see below.

struct HeightsProblem{TY <: AbstractVector, Tν <: Real}
    "Observations."
    y::TY
    "Degrees of freedom for prior on sigma."
    ν::Tν
end;

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

function (problem::HeightsProblem)(θ)
    @unpack y, ν = problem   # extract the data
    @unpack μ, σ = θ
    loglikelihood(Normal(μ, σ), y) + logpdf(TDist(ν), σ)
end;

Setup problem with data and inits.

obs = convert(Vector{Float64}, df2[:height]);
p = HeightsProblem(obs, 1.0);
p((μ = 178, σ = 5.0,))

-5170.976519811121

Write a function to return properly dimensioned transformation.

problem_transformation(p::HeightsProblem) =
    as((σ = asℝ₊, μ  = as(Real, 100, 250)), )

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.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Int64}}}},Main.ex-m4.1d.HeightsProblem{Array{Float64,1},Float64}},ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:σ, :μ),Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Int64}}}},Main.ex-m4.1d.HeightsProblem{Array{Float64,1},Float64}}},Float64},Float64,2,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(LogDensityProblems, Symbol("##1#2")){TransformedLogDensity{TransformVariables.TransformTuple{NamedTuple{(:σ, :μ),Tuple{TransformVariables.ShiftedExp{true,Float64},TransformVariables.ScaledShiftedLogistic{Int64}}}},Main.ex-m4.1d.HeightsProblem{Array{Float64,1},Float64}}},Float64},Float64,2},1}}}}(ForwardDiff AD wrapper for TransformedLogDensity of dimension 2, w/ chunk size 2)

Tune and sample.

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

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([2.0286197159657395, -0.5598009580347787], -1219.3532902020388, 2, DynamicHMC.AdjacentTurn, 0.999355966571135, 7), NUTS_Transition{Array{Float64,1},Float64}([2.1035408803379507, -0.5449044291088926], -1220.8308867900455, 2, DynamicHMC.DoubledTurn, 0.8087931551449169, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0047960713908113, -0.5834140026520368], -1224.7478851842423, 2, DynamicHMC.DoubledTurn, 0.6690022450147763, 3), NUTS_Transition{Array{Float64,1},Float64}([2.021966871204559, -0.5755147679551557], -1221.7039764132219, 1, DynamicHMC.DoubledTurn, 1.0, 1), NUTS_Transition{Array{Float64,1},Float64}([2.0681216020953186, -0.5411973943155898], -1220.3905216457963, 2, DynamicHMC.AdjacentTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([2.023825968604533, -0.5627733294098167], -1221.5327717007956, 1, DynamicHMC.AdjacentTurn, 0.8975386735329957, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0791311884297796, -0.5535537435369932], -1222.500391146448, 2, DynamicHMC.DoubledTurn, 0.6580120040050063, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0084870670329393, -0.5497362702975317], -1220.7896701433892, 2, DynamicHMC.DoubledTurn, 0.8947582098692664, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0103869662508096, -0.5723458578756284], -1221.5534073332535, 2, DynamicHMC.DoubledTurn, 0.8495108377157523, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0053232666801177, -0.5680308759022042], -1220.3760563574738, 2, DynamicHMC.DoubledTurn, 1.0, 3)  …  NUTS_Transition{Array{Float64,1},Float64}([2.085191497359603, -0.5652933004750212], -1223.3218878558084, 2, DynamicHMC.DoubledTurn, 0.5404595971858003, 3), NUTS_Transition{Array{Float64,1},Float64}([2.007096840858527, -0.549225247385076], -1220.0030560118958, 2, DynamicHMC.DoubledTurn, 0.9949901859986762, 3), NUTS_Transition{Array{Float64,1},Float64}([1.9746361949615314, -0.5604937238569611], -1221.360203163533, 2, DynamicHMC.DoubledTurn, 0.9234860301091139, 3), NUTS_Transition{Array{Float64,1},Float64}([2.048424483073958, -0.5476122787239495], -1221.3937069240772, 2, DynamicHMC.DoubledTurn, 0.954937282955858, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0479100635805856, -0.5648214083080484], -1220.398579456042, 2, DynamicHMC.DoubledTurn, 0.8972472907730523, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0387225892863428, -0.5560451762358397], -1219.7789856286995, 2, DynamicHMC.DoubledTurn, 0.936408840825818, 3), NUTS_Transition{Array{Float64,1},Float64}([1.9513439265413752, -0.5573030903257915], -1222.8752213941646, 1, DynamicHMC.AdjacentTurn, 0.4780338477943813, 3), NUTS_Transition{Array{Float64,1},Float64}([2.0726071337206804, -0.5477283134366904], -1224.3383152754377, 1, DynamicHMC.AdjacentTurn, 0.746636698100699, 3), NUTS_Transition{Array{Float64,1},Float64}([2.044106939886548, -0.5407054018673291], -1220.9618601719837, 2, DynamicHMC.DoubledTurn, 0.9022056832395459, 3), NUTS_Transition{Array{Float64,1},Float64}([2.060939497796793, -0.5795610748375615], -1222.5706438155003, 2, DynamicHMC.DoubledTurn, 0.7191095819843493, 3)], NUTS sampler in 2 dimensions
  stepsize (ϵ) ≈ 0.967
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.03805148661809207, 0.011042191267764418]
)

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

posterior = TransformVariables.transform.(Ref(problem_transformation(p)), get_position.(chain));

1000-element Array{NamedTuple{(:σ, :μ),Tuple{Float64,Float64}},1}:
 (σ = 7.6035840067253115, μ = 154.539027315885) 
 (σ = 8.195136596563385, μ = 155.05711396182065)
 (σ = 7.424579658161287, μ = 153.7221115915155) 
 (σ = 7.553166438784449, μ = 153.99479457711467)
 (σ = 7.9099511188486655, μ = 155.18636227271)  
 (σ = 7.567221571754781, μ = 154.43590119412184)
 (σ = 7.997517559236255, μ = 154.75604605444917)
 (σ = 7.45203438579522, μ = 154.88883976902596) 
 (σ = 7.466205958113087, μ = 154.1043560171591) 
 (σ = 7.428494893538293, μ = 154.25369709983994)
 ⋮                                              
 (σ = 7.441681570584242, μ = 154.906626428534)  
 (σ = 7.203998326135102, μ = 154.5149843864077) 
 (σ = 7.75567228252298, μ = 154.9627833689666)  
 (σ = 7.751683639522457, μ = 154.3648921716421) 
 (σ = 7.680791404306823, μ = 154.66945337686164)
 (σ = 7.038139970616106, μ = 154.62575528815503)
 (σ = 7.9455111479426, μ = 154.95874270842498)  
 (σ = 7.722259016011593, μ = 155.2035254611476) 
 (σ = 7.853344545400182, μ = 153.8550394887367) 

Extract the parameter posterior means: β,

posterior_μ = mean(last, posterior)

154.59032105775808

then σ:

posterior_σ = mean(first, posterior)

7.748504398651419

Effective sample sizes (of untransformed draws)

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

1×2 Array{Float64,2}:
 831.564  981.145

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
sigma   7.7641872 0.29928194 0.004732063 0.0055677898 1000
   mu 154.6055177 0.41989355 0.006639100 0.0085038356 1000

Quantiles:
         2.5%      25.0%       50.0%      75.0%       97.5%
sigma   7.21853   7.5560625   7.751355   7.9566775   8.410391
   mu 153.77992 154.3157500 154.602000 154.8820000 155.431000
";

"\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\nsigma   7.7641872 0.29928194 0.004732063 0.0055677898 1000\n   mu 154.6055177 0.41989355 0.006639100 0.0085038356 1000\n\nQuantiles:\n         2.5%      25.0%       50.0%      75.0%       97.5%\nsigma   7.21853   7.5560625   7.751355   7.9566775   8.410391\n   mu 153.77992 154.3157500 154.602000 154.8820000 155.431000\n"

Extract the parameter posterior means: β,

[posterior_μ, posterior_σ]

2-element Array{Float64,1}:
 154.59032105775808 
   7.748504398651419

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