Estimate Binomial draw probabilility

using DynamicHMCModels

Define a structure to hold the data.

struct BernoulliProblem
    "Total number of draws in the data."
    n::Int
    "Number of draws `==1` in the data"
    s::Vector{Int}
end;

Make the type callable with the parameters as a single argument.

function (problem::BernoulliProblem)((α, )::NamedTuple{(:α, )})
    @unpack n, s = problem        # extract the data
    loglikelihood(Binomial(n, α), s)
end

Create the data and complete setting up the problem.

obs = rand(Binomial(9, 2/3), 1)
p = BernoulliProblem(9, obs)
p((α = 0.5, ))

-1.4020427180880297

Write a function to return properly dimensioned transformation.

problem_transformation(p::BernoulliProblem) =
    as((α = as𝕀, ),  )

problem_transformation (generic function with 1 method)

Use a flat priors (the default, omitted) for α

P = TransformedLogDensity(problem_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P);

#import Zygote
#∇P = ADgradient(:Zygote, P);

Sample

chain, NUTS_tuned = NUTS_init_tune_mcmc(∇P, 1000)

(NUTS_Transition{Array{Float64,1},Float64}[NUTS_Transition{Array{Float64,1},Float64}([-1.016922479913577], -3.6541274324595574, 2, DynamicHMC.AdjacentTurn, 0.9124966260039936, 7), NUTS_Transition{Array{Float64,1},Float64}([-1.1645638780220735], -4.059328872229385, 1, DynamicHMC.DoubledTurn, 0.9466891548535735, 1), NUTS_Transition{Array{Float64,1},Float64}([0.15940491125123024], -3.758395824570598, 2, DynamicHMC.AdjacentTurn, 1.0, 7), NUTS_Transition{Array{Float64,1},Float64}([-0.3576996218592686], -3.3541649131014077, 2, DynamicHMC.DoubledTurn, 0.9461739470411451, 3), NUTS_Transition{Array{Float64,1},Float64}([0.5444872285799939], -4.005115368887483, 1, DynamicHMC.AdjacentTurn, 0.8154464316916088, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.30135496918986016], -3.7486528761570113, 1, DynamicHMC.AdjacentTurn, 0.962415822425741, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.5404634180568781], -2.981801219918469, 2, DynamicHMC.DoubledTurn, 0.9745912674970318, 3), NUTS_Transition{Array{Float64,1},Float64}([0.07563336698986356], -3.2295665159818667, 1, DynamicHMC.AdjacentTurn, 0.9676914252790375, 3), NUTS_Transition{Array{Float64,1},Float64}([-1.0946238294401176], -3.812365020723312, 2, DynamicHMC.DoubledTurn, 0.8861414854003081, 3), NUTS_Transition{Array{Float64,1},Float64}([-1.0246405874094588], -3.9168628313708553, 1, DynamicHMC.DoubledTurn, 1.0, 1)  …  NUTS_Transition{Array{Float64,1},Float64}([-0.3822466568968601], -2.811852938305918, 1, DynamicHMC.AdjacentTurn, 0.9881013676229246, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.14477322951527954], -2.8211756529395973, 1, DynamicHMC.AdjacentTurn, 0.9966308713812483, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.08138014620247075], -3.1474509395494947, 2, DynamicHMC.DoubledTurn, 0.9428735439491662, 3), NUTS_Transition{Array{Float64,1},Float64}([-1.0841070992761512], -3.8150943647166087, 2, DynamicHMC.DoubledTurn, 0.8719370768920228, 3), NUTS_Transition{Array{Float64,1},Float64}([0.24528450295798265], -5.03355389015074, 1, DynamicHMC.AdjacentTurn, 0.8970438564320938, 3), NUTS_Transition{Array{Float64,1},Float64}([0.47745328256923114], -3.3420991161287077, 2, DynamicHMC.DoubledTurn, 0.9763685922138214, 3), NUTS_Transition{Array{Float64,1},Float64}([0.47745328256923114], -4.01444021266594, 1, DynamicHMC.DoubledTurn, 0.8454214415932366, 1), NUTS_Transition{Array{Float64,1},Float64}([0.4865087898585469], -3.4770453136027366, 1, DynamicHMC.DoubledTurn, 0.9966223264496923, 1), NUTS_Transition{Array{Float64,1},Float64}([-1.2366661395167877], -6.87402131162114, 1, DynamicHMC.AdjacentTurn, 0.5925825853237678, 3), NUTS_Transition{Array{Float64,1},Float64}([-0.943951051310396], -4.060580888336315, 1, DynamicHMC.DoubledTurn, 1.0, 1)], NUTS sampler in 1 dimensions
  stepsize (ϵ) ≈ 0.933
  maximum depth = 10
  Gaussian kinetic energy, √diag(M⁻¹): [0.6060977800908485]
)

To get the posterior for $α$ use get_position and then transform back.

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

1000-element Array{NamedTuple{(:α,),Tuple{Float64}},1}:
 (α = 0.26562729747081265,)
 (α = 0.23783899406558046,)
 (α = 0.5397670569556116,) 
 (α = 0.41151653645788216,)
 (α = 0.632855642678423,)  
 (α = 0.4252262831657244,) 
 (α = 0.3680797859995336,) 
 (α = 0.5188993332843538,) 
 (α = 0.25074858153462365,)
 (α = 0.26412445656198774,)
 ⋮                         
 (α = 0.4638697759024505,) 
 (α = 0.4796661843220961,) 
 (α = 0.25272957339729346,)
 (α = 0.5610155181058538,) 
 (α = 0.617146324471203,)  
 (α = 0.617146324471203,)  
 (α = 0.6192836482804898,) 
 (α = 0.22501682547729007,)
 (α = 0.2801029399583933,) 

Extract the parameter.

posterior_α = first.(posterior);

1000-element Array{Float64,1}:
 0.26562729747081265
 0.23783899406558046
 0.5397670569556116 
 0.41151653645788216
 0.632855642678423  
 0.4252262831657244 
 0.3680797859995336 
 0.5188993332843538 
 0.25074858153462365
 0.26412445656198774
 ⋮                  
 0.4638697759024505 
 0.4796661843220961 
 0.25272957339729346
 0.5610155181058538 
 0.617146324471203  
 0.617146324471203  
 0.6192836482804898 
 0.22501682547729007
 0.2801029399583933 

check the effective sample size

ess_α = effective_sample_size(posterior_α)

352.0211439484193

NUTS-specific statistics

NUTS_statistics(chain)

Hamiltonian Monte Carlo sample of length 1000
  acceptance rate mean: 0.94, min/25%/median/75%/max: 0.28 0.92 0.98 1.0 1.0
  termination: AdjacentTurn => 42% DoubledTurn => 58%
  depth: 1 => 64% 2 => 35% 3 => 0%

check the mean

mean(posterior_α)

0.4565266715691886

This page was generated using Literate.jl.