Inference example for VFTS 243
In this example we try to match the observed properties of VFTS 243. We start by loading up the SideKicks package, as well the Distributions package which we use to define the priors.
using SideKicks
using DistributionsWe set up our observations using the @Observations macro. This macro will return a tuple containing an instance of Observation, as well as a string with the code needed to initialize the Observation instance. The reason for this is that when we save results, we also want to save the precise observations that were used for the MCMC.
Currently, The likelihoods are taken to be Normal distributions. This can potentially be generalized in the future to arbitrary user-defined likelihood functions.
obs = SideKicks.@Observations([
[:P_f, 10.4031, 0.01, day],
[:e_f, 0.017, 0.012, 1],
[:m1_f, 25.0, 2.3, m_sun],
[:K1, 81.4, 1.3, km_per_s],
[:v_N, 138.8, 7.6, km_per_s], # Gaia 4' w/ MCMC
[:v_E, 409.3, 9.3, km_per_s], # Gaia 4' w/ MCMC
[:v_r, 261.5, 0.42, km_per_s], # Almeida
[:ω_f, 66, 53, degree]
])(SideKicks.Observations([:P_f, :e_f, :m1_f, :K1, :v_N, :v_E, :v_r, :ω_f], [10.4031, 0.017, 25.0, 81.4, 138.8, 409.3, 261.5, 66.0], [0.01, 0.012, 2.3, 1.3, 7.6, 9.3, 0.42, 53.0], [86400.0, 1.0, 1.9884098706980504e33, 100000.0, 100000.0, 100000.0, 100000.0, 0.017453292519943295]), "(SideKicks.Observations)([[:P_f, 10.4031, 0.01, day], [:e_f, 0.017, 0.012, 1], [:m1_f, 25.0, 2.3, m_sun], [:K1, 81.4, 1.3, km_per_s], [:v_N, 138.8, 7.6, km_per_s], [:v_E, 409.3, 9.3, km_per_s], [:v_r, 261.5, 0.42, km_per_s], [:ω_f, 66, 53, degree]])")Next, we choose our priors. In here we are using distributions defined in the Distributions.jl` package.
priors = SideKicks.@Priors(
logm1_dist = Uniform(0.1,3), # in log(Msun)
logm2_dist = Uniform(0.1,3), # in log(Msun)
logP_dist = Uniform(-1,3), # in log(days)
vkick_dist = Uniform(0,4), # in 100 km/s
frac_dist = Uniform(0,1.0),
e_dist = Uniform(0,0.01),
venv_N_100kms_dist = Normal(146/100, 40/100), # Gaia 4' w/ MCMC
venv_E_100kms_dist = Normal(393/100, 42/100), # Gaia 4' w/ MCMC
venv_r_100kms_dist = Normal(270.3/100, 11.1/100) # Almeida w/ MCMC
)(SideKicks.Priors(Distributions.Uniform{Float64}(a=0.1, b=3.0), Distributions.Uniform{Float64}(a=0.1, b=3.0), Distributions.Uniform{Float64}(a=-1.0, b=3.0), Distributions.Uniform{Float64}(a=0.0, b=4.0), Distributions.Uniform{Float64}(a=0.0, b=1.0), Distributions.Uniform{Float64}(a=0.0, b=0.01), Distributions.Normal{Float64}(μ=1.46, σ=0.4), Distributions.Normal{Float64}(μ=3.93, σ=0.42), Distributions.Normal{Float64}(μ=2.7030000000000003, σ=0.111), missing, missing, missing, missing), "SideKicks.Priors(logm1_dist = Uniform(0.1, 3), logm2_dist = Uniform(0.1, 3), logP_dist = Uniform(-1, 3), vkick_dist = Uniform(0, 4), frac_dist = Uniform(0, 1.0), e_dist = Uniform(0, 0.01), venv_N_100kms_dist = Normal(146 / 100, 40 / 100), venv_E_100kms_dist = Normal(393 / 100, 42 / 100), venv_r_100kms_dist = Normal(270.3 / 100, 11.1 / 100))")Finally, we run the MCMC using the NUTS algorithm. In the example below we are computing 3 chains with 100 samples each. This in practice is very low, and a ballpark suggestion would be to total at least $\sim10^4$ samples for any meaningful analysis, and an order of magnitude higher for final science runs.
kick_mcmc = SideKicks.KickMCMC(
which_model = :general,
observations = obs,
priors = priors,
nuts_warmup_count = 500,
nuts_acceptance_rate = 0.8,
nsamples = 1000,
nchains = 4)SideKicks.KickMCMC(DynamicPPL.Model{SideKicks.var"#create_mcmc_model#15"{Symbol, typeof(SideKicks.arbitraryEjectaBH), Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}}, (:props, :obs_vals, :obs_errs), (), (), Tuple{Vector{Symbol}, Vector{Float64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}(SideKicks.var"#create_mcmc_model#15"{Symbol, typeof(SideKicks.arbitraryEjectaBH), Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}}(:Cauchy, SideKicks.arbitraryEjectaBH, Distributions.Normal{Float64}(μ=2.7030000000000003, σ=0.111), Distributions.Normal{Float64}(μ=3.93, σ=0.42), Distributions.Normal{Float64}(μ=1.46, σ=0.4), Distributions.Uniform{Float64}(a=0.0, b=1.0), Distributions.Uniform{Float64}(a=0.0, b=4.0), Distributions.Uniform{Float64}(a=0.0, b=0.01), Distributions.Uniform{Float64}(a=-1.0, b=3.0), Distributions.Uniform{Float64}(a=0.1, b=3.0), Distributions.Uniform{Float64}(a=0.1, b=3.0), Core.Box(SideKicks.var"#create_mcmc_model#15"{Symbol, typeof(SideKicks.arbitraryEjectaBH), Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}}(#= circular reference @-2 =#))), (props = [:P_f, :e_f, :m1_f, :K1, :v_N, :v_E, :v_r, :ω_f], obs_vals = [898827.84, 0.017, 4.971024676745126e34, 8.140000000000001e6, 1.3880000000000002e7, 4.093e7, 2.615e7, 1.1519173063162575], obs_errs = [864.0, 0.012, 4.5733427026055155e33, 130000.0, 760000.0, 930000.0000000001, 42000.0, 0.9250245035569946]), NamedTuple(), DynamicPPL.DefaultContext()), DynamicPPL.Model{SideKicks.var"#create_mcmc_model#15"{Symbol, typeof(SideKicks.arbitraryEjectaBH), Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}}, (:props, :obs_vals, :obs_errs), (), (), Tuple{Vector{Symbol}, Vector{Float64}, Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}(SideKicks.var"#create_mcmc_model#15"{Symbol, typeof(SideKicks.arbitraryEjectaBH), Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}}(:Normal, SideKicks.arbitraryEjectaBH, Distributions.Normal{Float64}(μ=2.7030000000000003, σ=0.111), Distributions.Normal{Float64}(μ=3.93, σ=0.42), Distributions.Normal{Float64}(μ=1.46, σ=0.4), Distributions.Uniform{Float64}(a=0.0, b=1.0), Distributions.Uniform{Float64}(a=0.0, b=4.0), Distributions.Uniform{Float64}(a=0.0, b=0.01), Distributions.Uniform{Float64}(a=-1.0, b=3.0), Distributions.Uniform{Float64}(a=0.1, b=3.0), Distributions.Uniform{Float64}(a=0.1, b=3.0), Core.Box(SideKicks.var"#create_mcmc_model#15"{Symbol, typeof(SideKicks.arbitraryEjectaBH), Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Normal{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}, Distributions.Uniform{Float64}}(#= circular reference @-2 =#))), (props = [:P_f, :e_f, :m1_f, :K1, :v_N, :v_E, :v_r, :ω_f], obs_vals = [898827.84, 0.017, 4.971024676745126e34, 8.140000000000001e6, 1.3880000000000002e7, 4.093e7, 2.615e7, 1.1519173063162575], obs_errs = [864.0, 0.012, 4.5733427026055155e33, 130000.0, 760000.0, 930000.0000000001, 42000.0, 0.9250245035569946]), NamedTuple(), DynamicPPL.DefaultContext()), SideKicks.Observations([:P_f, :e_f, :m1_f, :K1, :v_N, :v_E, :v_r, :ω_f], [10.4031, 0.017, 25.0, 81.4, 138.8, 409.3, 261.5, 66.0], [0.01, 0.012, 2.3, 1.3, 7.6, 9.3, 0.42, 53.0], [86400.0, 1.0, 1.9884098706980504e33, 100000.0, 100000.0, 100000.0, 100000.0, 0.017453292519943295]), "(SideKicks.Observations)([[:P_f, 10.4031, 0.01, day], [:e_f, 0.017, 0.012, 1], [:m1_f, 25.0, 2.3, m_sun], [:K1, 81.4, 1.3, km_per_s], [:v_N, 138.8, 7.6, km_per_s], [:v_E, 409.3, 9.3, km_per_s], [:v_r, 261.5, 0.42, km_per_s], [:ω_f, 66, 53, degree]])", SideKicks.Priors(Distributions.Uniform{Float64}(a=0.1, b=3.0), Distributions.Uniform{Float64}(a=0.1, b=3.0), Distributions.Uniform{Float64}(a=-1.0, b=3.0), Distributions.Uniform{Float64}(a=0.0, b=4.0), Distributions.Uniform{Float64}(a=0.0, b=1.0), Distributions.Uniform{Float64}(a=0.0, b=0.01), Distributions.Normal{Float64}(μ=1.46, σ=0.4), Distributions.Normal{Float64}(μ=3.93, σ=0.42), Distributions.Normal{Float64}(μ=2.7030000000000003, σ=0.111), missing, missing, missing, missing), "SideKicks.Priors(logm1_dist = Uniform(0.1, 3), logm2_dist = Uniform(0.1, 3), logP_dist = Uniform(-1, 3), vkick_dist = Uniform(0, 4), frac_dist = Uniform(0, 1.0), e_dist = Uniform(0, 0.01), venv_N_100kms_dist = Normal(146 / 100, 40 / 100), venv_E_100kms_dist = Normal(393 / 100, 42 / 100), venv_r_100kms_dist = Normal(270.3 / 100, 11.1 / 100))", Dict(:P_f => [899631.9031002622 899128.2365316936 899567.7176608172 899274.2186061341; 897420.9590948492 898364.0776499602 899280.8119107054 898661.5158801406; … ; 899853.3482188154 902558.2047020107 898875.7129623046 896802.1472756597; 899667.0430973007 897347.7153521972 898940.5618871566 907406.0788773461], :venv_N => [1.3876483202127969e7 1.2621530687144091e7 1.4068447071361382e7 1.2572763598337747e7; 1.4297599649547035e7 1.2681803314372845e7 1.4013590948218038e7 1.5102794465687463e7; … ; 1.5956360603457477e7 1.7627150018358894e7 1.2305354307782372e7 9.722302898470417e6; 1.4216782827237887e7 1.7964819484936602e7 1.393845812017831e7 9.90533834806879e6], :ϕ => [1.1855987607820342 0.34423732021334685 0.22719236136917054 6.059693841940314; 1.1800452040065998 0.2710937181347743 4.247340964941075 4.9401887733347305; … ; 2.4468645221823544 2.295664622174951 1.2973145884130466 0.6384852168876669; 2.107615085076962 2.131552186127755 1.6707341096559 0.6983512516108324], :vsys => [2.4459897203830904e6 4.74591727359189e6 810096.3371508452 164273.80881000054; 2.252121160849232e6 4.598989648369307e6 466188.0807210991 1.0971275371248862e6; … ; 2.1120811611542227e6 3.5189278072270844e6 2.4411074574470157e6 1.7674036107103191e6; 1.7807666745114257e6 3.206321547812187e6 1.4435516203253793e6 1.8567173424360012e6], :a_f => [3.974197911516519e12 4.3219480324469956e12 5.336350505940822e12 4.187743135025976e12; 4.07599781431358e12 4.2924732257835117e12 4.449246065584152e12 4.2589297315169077e12; … ; 3.9624942503117715e12 4.527507089729218e12 4.925581015847777e12 4.752289113594853e12; 4.418765925615615e12 4.49279148936757e12 3.917708057959593e12 4.866225974013546e12], :m2_i => [4.5420429663070215e33 1.0255206550372551e34 6.529789111609211e34 9.055041353242428e33; 4.472022577130002e33 9.72492981954303e33 1.8589692075620655e34 4.613751095292726e33; … ; 5.375767771195607e33 2.112189025645653e34 2.853017619425047e34 3.222067675326792e34; 5.501945541574509e33 1.9356664733722574e34 1.6782648171627284e34 3.581508898025151e34], :m2_f => [3.658346955348851e32 1.0189180748545276e34 6.006584974215178e34 8.982047826105577e33; 3.4680697658736314e32 9.639946375569758e33 1.7325265225037192e34 2.4240092352998267e33; … ; 1.6777356685227182e33 1.7502680402613544e34 2.845474641511032e34 2.685617910056339e34; 1.664468978164774e33 1.6287018641077607e34 1.5345771063226315e34 2.9765933782810635e34], :e_i => [0.003283606870221012 0.0015448147167173942 0.00238191766586427 0.0019022910718393565; 0.003370977828596793 0.0017171465728704183 0.007984786172564953 0.00876975643743151; … ; 0.0079152573677324 0.006566989623737115 0.007984928979181561 0.001226729868378267; 0.007824588410649536 0.006410373704541231 0.0026614356894650287 0.00279881089829437], :vsys_r => [-1.114922850497939e6 -1.2981358022188311e6 -62372.49692616431 -10262.22165019547; -805673.2861609499 -1.6713906647240615e6 -79873.29706496377 -367949.3551159365; … ; -532037.7174065365 -1.8906511514613263e6 -1.6932254290861562e6 -913211.4145697725; -421278.13856933184 -1.521111855904251e6 -861526.348945061 -964913.2511110031], :K2 => [1.830420922286099e7 3.800748217623164e7 7.439969989771869e6 2.3018731464705925e7; 2.1233625634853058e7 3.7718798881862886e7 2.223512577914477e7 2.0512358588310216e7; … ; 9.9445515838298e6 2.4050177549604084e7 1.7427862623049874e7 1.5631269031815333e7; 1.1236761356232785e7 2.419777516199235e7 1.5328656177662924e7 1.4676130614657408e7]…), MCMC chain (1000×31×4 Array{Float64, 3}), 500, 0.8, 1000, Dict(:P_f => 4176.165226362999, :venv_N => 138.28239556854203, :v_r => 3518.0964044340976, :vsys => 17.3480710278482, :a_f => 28.824220868751226, :m2_i => 13.834183346058104, :m2_f => 12.185202576460457, :e_i => 37.5077489258474, :vsys_r => 27.872405218894887, :K2 => 14.293181560986527…), Dict(:P_f => 1.0035902366146248, :venv_N => 1.0306699918103055, :v_r => 1.0083754407418068, :vsys => 1.340889606629902, :a_f => 1.2019148561293034, :m2_i => 1.542710465554965, :m2_f => 1.7386735129863196, :e_i => 1.123414614732654, :vsys_r => 1.1459085197754773, :K2 => 1.4875062807392396…))Results from the MCMC can be saved to an HDF5 file. In practice, it is ideal to separate the computation of the MCMC and its analysis, so plotting of the results is done separetely.
SideKicks.SaveResults("vfts243_results.hdf5", kick_mcmc)This page was generated using Literate.jl.