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);| 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);| 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 | 161.29 | 48.9879 | 39.0 | 1 |
| 20 | 156.21 | 42.7227 | 29.0 | 0 |
| 21 | 146.4 | 35.4936 | 56.0 | 1 |
| 22 | 148.59 | 37.9033 | 45.0 | 0 |
| 23 | 147.32 | 35.4652 | 19.0 | 0 |
| 24 | 147.955 | 40.313 | 29.0 | 1 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
Show the first six rows of the dataset.
first(df2, 6)| 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 |
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.102069832888Write a function to return properly dimensioned transformation.
problem_transformation(p::ConstraintHeightsProblem) =
as((μ = as(Real, 100, 250), σ = asℝ₊), )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.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.TransformNamedTuple{(:μ, :σ),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.TransformNamedTuple{(:μ, :σ),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);MCMC, adapting ϵ (75 steps)
0.00044 s/step ...done
MCMC, adapting ϵ (25 steps)
0.00027 s/step ...done
MCMC, adapting ϵ (50 steps)
0.00037 s/step ...done
MCMC, adapting ϵ (100 steps)
7.0e-5 s/step ...done
MCMC, adapting ϵ (200 steps)
0.00019 s/step ...done
MCMC, adapting ϵ (400 steps)
5.2e-5 s/step ...done
MCMC, adapting ϵ (50 steps)
6.6e-5 s/step ...done
MCMC (1000 steps)
5.2e-5 s/step ...done
(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([0.07643999083058024, 3.232311770507498], -1623.5409215634847, 1, DynamicHMC.AdjacentTurn, 0.9327842961222735, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07559166423535812, 3.2047259650136906], -1622.657695437832, 2, DynamicHMC.AdjacentTurn, 0.997095193602199, 5), NUTS_Transition{Array{Float64,1},Float64}([0.0743870151663278, 3.1739506828115833], -1623.3144625273417, 2, DynamicHMC.AdjacentTurn, 0.9236880790861411, 5), NUTS_Transition{Array{Float64,1},Float64}([0.0754380737781423, 3.1710479299878713], -1622.8125575723511, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07560069385136872, 3.234600571767246], -1623.123356239427, 3, DynamicHMC.DoubledTurn, 0.9865438383597953, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07559754388445492, 3.2252414427501566], -1622.8976135198027, 2, DynamicHMC.DoubledTurn, 0.9898530393346631, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07188387995012636, 3.172854292641694], -1625.040850635383, 3, DynamicHMC.DoubledTurn, 0.8083146699630536, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07229059864128068, 3.150430871680783], -1624.7128456570806, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07138269255688302, 3.1517611346804846], -1624.904603604297, 1, DynamicHMC.DoubledTurn, 0.8934776864843658, 1), NUTS_Transition{Array{Float64,1},Float64}([0.08098992050713445, 3.164749985449015], -1625.5087191167393, 2, DynamicHMC.DoubledTurn, 0.9895443528434789, 3) … NUTS_Transition{Array{Float64,1},Float64}([0.0772613775877543, 3.1866533825217265], -1622.8493672360112, 3, DynamicHMC.DoubledTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07487838873055738, 3.21906859969851], -1622.9690815103093, 2, DynamicHMC.DoubledTurn, 0.9427987351762085, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07310859786668705, 3.2070292395343976], -1623.23372676545, 2, DynamicHMC.DoubledTurn, 0.9554631470558702, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07138258756014565, 3.175958474625386], -1624.749880295624, 2, DynamicHMC.AdjacentTurn, 0.9440174545046279, 7), NUTS_Transition{Array{Float64,1},Float64}([0.07670991137217636, 3.232912352874054], -1624.446325463207, 2, DynamicHMC.AdjacentTurn, 0.9934779676997267, 7), NUTS_Transition{Array{Float64,1},Float64}([0.08046315855045277, 3.2483868154696807], -1624.3823451493558, 2, DynamicHMC.DoubledTurn, 0.7883629247830379, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07960243245813693, 3.239175900466854], -1623.971205961067, 1, DynamicHMC.DoubledTurn, 1.0, 1), NUTS_Transition{Array{Float64,1},Float64}([0.07219887445831179, 3.1576186886162336], -1624.8104251251625, 2, DynamicHMC.DoubledTurn, 0.898195332828637, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07253870694120161, 3.1461184068370227], -1624.5187920493095, 2, DynamicHMC.DoubledTurn, 1.0, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07981584807833692, 3.167943742023101], -1625.4827456193746, 2, DynamicHMC.DoubledTurn, 0.9949570589444431, 3)], NUTS sampler in 2 dimensions
stepsize (ϵ) ≈ 0.651
maximum depth = 10
Gaussian kinetic energy, √diag(M⁻¹): [0.0037964059009868573, 0.03699532022253563]
)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.86510470751574, σ = 25.3381653392911)
(μ = 177.83333837246704, σ = 24.648744472516753)
(μ = 177.7882274826198, σ = 23.90172621292757)
(μ = 177.8275869336099, σ = 23.83244600996988)
(μ = 177.83367649975636, σ = 25.3962257831015)
(μ = 177.8335585446125, σ = 25.15964803715661)
(μ = 177.69448533195552, σ = 23.8755349558064)
(μ = 177.7097174855122, σ = 23.34612159679764)
(μ = 177.6757148944217, σ = 23.37719874435131)
(μ = 178.03546297401874, σ = 23.682822242426052)
⋮
(μ = 177.80662835588933, σ = 25.004819835314272)
(μ = 177.7403519606342, σ = 24.70558272955475)
(μ = 177.67571096205552, σ = 23.949764112066322)
(μ = 177.87521190411826, σ = 25.35338756523745)
(μ = 178.0157415482289, σ = 25.748768879699515)
(μ = 177.98351595105134, σ = 25.51268808728517)
(μ = 177.70628231290212, σ = 23.51453377787187)
(μ = 177.71900936005, σ = 23.245659044480533)
(μ = 177.9915063382548, σ = 23.758580323695888)Extract the parameter posterior means: μ,
posterior_μ = mean(first, posterior)177.86595784875962Extract the parameter posterior means: μ,
posterior_σ = mean(last, posterior)24.56977882979769Effective sample sizes (of untransformed draws)
ess = mapslices(effective_sample_size,
get_position_matrix(chain); dims = 1)NUTS-specific statistics
NUTS_statistics(chain)Hamiltonian Monte Carlo sample of length 1000
acceptance rate mean: 0.93, min/25%/median/75%/max: 0.38 0.9 0.96 1.0 1.0
termination: AdjacentTurn => 37% DoubledTurn => 63%
depth: 1 => 17% 2 => 68% 3 => 15% 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.86595784875962
24.56977882979769end of m4.5d.jl#- This page was generated using Literate.jl.