Heights problem with restricted prior on mu.

Result is not conform cmdstan result

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

No covariates, just height observations.

struct ConstraintHeightsProblem{TY <: AbstractVector}
    "Observations."
    y::TY
end;

Very constraint prior on μ. Flat σ.

function (problem::ConstraintHeightsProblem)(θ)
    @unpack y = problem   # extract the data
    @unpack μ, σ = θ
    loglikelihood(Normal(μ, σ), y) + logpdf(Normal(178, 0.1), μ) +
    logpdf(Uniform(0, 50), σ)
end;

Define problem with data and inits.

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

-5169.102069832888

Write a function to return properly dimensioned transformation.

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

problem_transformation (generic function with 1 method)

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

FSample from the posterior.

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

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([0.07523901131124944, 3.161603349484391], -1624.2311157293077, 3, DynamicHMC.DoubledTurn, 0.9963483729570839, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07873799640185371, 3.2050223947873504], -1623.4016993246387, 2, DynamicHMC.AdjacentTurn, 0.9685691412718824, 5), NUTS_Transition{Array{Float64,1},Float64}([0.07873799640185371, 3.2050223947873504], -1625.084223476229, 2, DynamicHMC.DoubledTurn, 0.7383123895838916, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07801811574558898, 3.2035618476418892], -1622.6278822119937, 1, DynamicHMC.DoubledTurn, 1.0, 1), NUTS_Transition{Array{Float64,1},Float64}([0.07314503507186369, 3.2012317616803374], -1623.0425711911528, 2, DynamicHMC.DoubledTurn, 0.9531413216462439, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07928884534151467, 3.1449324919497896], -1624.2138010897954, 2, DynamicHMC.AdjacentTurn, 0.9790607712955194, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07478255986707702, 3.1546850213698137], -1624.1706542067236, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07328879171381861, 3.1638455061043906], -1624.2257184501113, 2, DynamicHMC.DoubledTurn, 0.8603324806073774, 3), NUTS_Transition{Array{Float64,1},Float64}([0.0794368405287834, 3.2140031363696093], -1623.5011811298912, 3, DynamicHMC.DoubledTurn, 0.9993090665045098, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07469924361534341, 3.2119271745459104], -1623.5975068889193, 2, DynamicHMC.DoubledTurn, 0.9289173332938367, 3)  …  NUTS_Transition{Array{Float64,1},Float64}([0.07905467133634535, 3.156536442845271], -1626.8595495418574, 2, DynamicHMC.AdjacentTurn, 0.6827757901214764, 5), NUTS_Transition{Array{Float64,1},Float64}([0.07808278690145559, 3.2026288467087944], -1624.386957882765, 2, DynamicHMC.AdjacentTurn, 0.9464858576311265, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07398550912164537, 3.1627958353118735], -1623.3230583877105, 2, DynamicHMC.AdjacentTurn, 0.9678116710700052, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07385305815673777, 3.14444723658436], -1627.2101978127894, 3, DynamicHMC.DoubledTurn, 0.8340339996920061, 7), NUTS_Transition{Array{Float64,1},Float64}([0.08164467987906972, 3.179935921242356], -1625.431953890554, 2, DynamicHMC.AdjacentTurn, 0.878907122293462, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07397465414730382, 3.1855905817733947], -1625.8464237752269, 1, DynamicHMC.AdjacentTurn, 0.8632194246219171, 3), NUTS_Transition{Array{Float64,1},Float64}([0.08386536788535154, 3.175586695774094], -1626.7856947446653, 2, DynamicHMC.DoubledTurn, 0.6835121707136295, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07445385378651319, 3.230395876078697], -1626.6619235526102, 2, DynamicHMC.AdjacentTurn, 0.9932587188797667, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07829492378219953, 3.1577993375457245], -1624.1744970185016, 2, DynamicHMC.AdjacentTurn, 0.935613857159054, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07412469085625858, 3.1864122350143944], -1623.3532930898539, 2, DynamicHMC.DoubledTurn, 1.0, 3)], NUTS sampler in 2 dimensions
  stepsize (ϵ) ≈ 0.574
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.004653982067163362, 0.034246180315165]
)

Undo 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}:
 (μ = 177.82013267351277, σ = 23.608418144910367)
 (μ = 177.95115034222383, σ = 24.65605217731887) 
 (μ = 177.95115034222383, σ = 24.65605217731887) 
 (μ = 177.92419623467694, σ = 24.62006713601456) 
 (μ = 177.74171653074487, σ = 24.562767046200854)
 (μ = 177.97177497007007, σ = 23.218108011175385)
 (μ = 177.80303979947413, σ = 23.445651049218046)
 (μ = 177.74710018652914, σ = 23.66141130322116) 
 (μ = 177.97731605992524, σ = 24.878479094734015)
 (μ = 177.79991979930142, σ = 24.82688589318895) 
 ⋮                                               
 (μ = 177.92661771331416, σ = 24.597107302842215)
 (μ = 177.77319170321886, σ = 23.636587641482024)
 (μ = 177.76823157073838, σ = 23.206844032861824)
 (μ = 178.0599759049541, σ = 24.04521271533949)  
 (μ = 177.7727851981407, σ = 24.181565381335744) 
 (μ = 178.14310928289905, σ = 23.9408617512333)  
 (μ = 177.79073046364852, σ = 25.289666563583772)
 (μ = 177.93456070029166, σ = 23.518782036936493)
 (μ = 177.77840386891347, σ = 24.2014424077971)  

Extract the parameter posterior means: μ,

posterior_μ = mean(first, posterior)

177.86194327932654

Extract the parameter posterior means: μ,

posterior_σ = mean(last, posterior)

24.50153416579062

Effective sample sizes (of untransformed draws)

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

1×2 Array{Float64,2}:
 1000.0  346.779

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.93, min/25%/median/75%/max: 0.3 0.9 0.96 0.99 1.0
  termination: AdjacentTurn => 40% DoubledTurn => 60%
  depth: 1 => 15% 2 => 68% 3 => 17% 4 => 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
sigma  24.604616 0.946911707 0.0149719887 0.0162406632 1000
   mu 177.864069 0.102284043 0.0016172527 0.0013514459 1000

Quantiles:
         2.5%       25.0%     50.0%     75.0%     97.5%
sigma  22.826377  23.942275  24.56935  25.2294  26.528368
   mu 177.665000 177.797000 177.86400 177.9310 178.066000
";

"\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  24.604616 0.946911707 0.0149719887 0.0162406632 1000\n   mu 177.864069 0.102284043 0.0016172527 0.0013514459 1000\n\nQuantiles:\n         2.5%       25.0%     50.0%     75.0%     97.5%\nsigma  22.826377  23.942275  24.56935  25.2294  26.528368\n   mu 177.665000 177.797000 177.86400 177.9310 178.066000\n"

Extract the parameter posterior means: β,

[posterior_μ, posterior_σ]

2-element Array{Float64,1}:
 177.86194327932654
  24.50153416579062

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