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);┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC ~/.julia/packages/AbstractMCMC/oqm6Y/src/sample.jl:544
Sampling (1 thread) 0%| | ETA: N/A
┌ Info: Found initial step size
└ ϵ = 0.4
Sampling (1 thread) 0%|▏ | ETA: 1:22:06
Sampling (1 thread) 1%|▎ | ETA: 0:41:13
Sampling (1 thread) 2%|▌ | ETA: 0:27:50
Sampling (1 thread) 2%|▋ | ETA: 0:21:09
Sampling (1 thread) 2%|▊ | ETA: 0:17:08
Sampling (1 thread) 3%|▉ | ETA: 0:14:26
Sampling (1 thread) 4%|█▏ | ETA: 0:12:31
Sampling (1 thread) 4%|█▎ | ETA: 0:11:05
Sampling (1 thread) 4%|█▍ | ETA: 0:09:57
Sampling (1 thread) 5%|█▌ | ETA: 0:09:03
Sampling (1 thread) 6%|█▊ | ETA: 0:08:19
Sampling (1 thread) 6%|█▉ | ETA: 0:07:42
Sampling (1 thread) 6%|██ | ETA: 0:07:10
Sampling (1 thread) 7%|██▏ | ETA: 0:06:42
Sampling (1 thread) 8%|██▍ | ETA: 0:06:19
Sampling (1 thread) 8%|██▌ | ETA: 0:05:58
Sampling (1 thread) 8%|██▋ | ETA: 0:05:42
Sampling (1 thread) 9%|██▊ | ETA: 0:05:26
Sampling (1 thread) 10%|███ | ETA: 0:05:11
Sampling (1 thread) 10%|███▏ | ETA: 0:04:58
Sampling (1 thread) 10%|███▎ | ETA: 0:04:47
Sampling (1 thread) 11%|███▍ | ETA: 0:04:36
Sampling (1 thread) 12%|███▋ | ETA: 0:04:26
Sampling (1 thread) 12%|███▊ | ETA: 0:04:17
Sampling (1 thread) 12%|███▉ | ETA: 0:04:08
Sampling (1 thread) 13%|████ | ETA: 0:04:00
Sampling (1 thread) 14%|████▏ | ETA: 0:03:53
Sampling (1 thread) 14%|████▍ | ETA: 0:03:46
Sampling (1 thread) 14%|████▌ | ETA: 0:03:40
Sampling (1 thread) 15%|████▋ | ETA: 0:03:34
Sampling (1 thread) 16%|████▊ | ETA: 0:03:28
Sampling (1 thread) 16%|█████ | ETA: 0:03:23
Sampling (1 thread) 16%|█████▏ | ETA: 0:03:17
Sampling (1 thread) 17%|█████▎ | ETA: 0:03:13
Sampling (1 thread) 18%|█████▍ | ETA: 0:03:08
Sampling (1 thread) 18%|█████▋ | ETA: 0:03:04
Sampling (1 thread) 18%|█████▊ | ETA: 0:03:00
Sampling (1 thread) 19%|█████▉ | ETA: 0:02:56
Sampling (1 thread) 20%|██████ | ETA: 0:02:52
Sampling (1 thread) 20%|██████▎ | ETA: 0:02:48
Sampling (1 thread) 20%|██████▍ | ETA: 0:02:45
Sampling (1 thread) 21%|██████▌ | ETA: 0:02:42
Sampling (1 thread) 22%|██████▋ | ETA: 0:02:39
Sampling (1 thread) 22%|██████▉ | ETA: 0:02:36
Sampling (1 thread) 22%|███████ | ETA: 0:02:33
Sampling (1 thread) 23%|███████▏ | ETA: 0:02:30
Sampling (1 thread) 24%|███████▎ | ETA: 0:02:27
Sampling (1 thread) 24%|███████▌ | ETA: 0:02:24
Sampling (1 thread) 24%|███████▋ | ETA: 0:02:22
Sampling (1 thread) 25%|███████▊ | ETA: 0:02:20
┌ Info: Found initial step size
└ ϵ = 0.22138671875
Sampling (1 thread) 26%|███████▉ | ETA: 0:02:46
Sampling (1 thread) 26%|████████ | ETA: 0:02:43
Sampling (1 thread) 26%|████████▎ | ETA: 0:02:40
Sampling (1 thread) 27%|████████▍ | ETA: 0:02:38
Sampling (1 thread) 28%|████████▌ | ETA: 0:02:35
Sampling (1 thread) 28%|████████▋ | ETA: 0:02:34
Sampling (1 thread) 28%|████████▉ | ETA: 0:02:32
Sampling (1 thread) 29%|█████████ | ETA: 0:02:29
Sampling (1 thread) 30%|█████████▏ | ETA: 0:02:27
Sampling (1 thread) 30%|█████████▎ | ETA: 0:02:25
Sampling (1 thread) 30%|█████████▌ | ETA: 0:02:22
Sampling (1 thread) 31%|█████████▋ | ETA: 0:02:20
Sampling (1 thread) 32%|█████████▊ | ETA: 0:02:18
Sampling (1 thread) 32%|█████████▉ | ETA: 0:02:16
Sampling (1 thread) 32%|██████████▏ | ETA: 0:02:13
Sampling (1 thread) 33%|██████████▎ | ETA: 0:02:11
Sampling (1 thread) 34%|██████████▍ | ETA: 0:02:09
Sampling (1 thread) 34%|██████████▌ | ETA: 0:02:07
Sampling (1 thread) 34%|██████████▊ | ETA: 0:02:05
Sampling (1 thread) 35%|██████████▉ | ETA: 0:02:04
Sampling (1 thread) 36%|███████████ | ETA: 0:02:02
Sampling (1 thread) 36%|███████████▏ | ETA: 0:02:00
Sampling (1 thread) 36%|███████████▍ | ETA: 0:01:58
Sampling (1 thread) 37%|███████████▌ | ETA: 0:01:56
Sampling (1 thread) 38%|███████████▋ | ETA: 0:01:55
Sampling (1 thread) 38%|███████████▊ | ETA: 0:01:53
Sampling (1 thread) 38%|███████████▉ | ETA: 0:01:51
Sampling (1 thread) 39%|████████████▏ | ETA: 0:01:50
Sampling (1 thread) 40%|████████████▎ | ETA: 0:01:48
Sampling (1 thread) 40%|████████████▍ | ETA: 0:01:47
Sampling (1 thread) 40%|████████████▌ | ETA: 0:01:45
Sampling (1 thread) 41%|████████████▊ | ETA: 0:01:44
Sampling (1 thread) 42%|████████████▉ | ETA: 0:01:42
Sampling (1 thread) 42%|█████████████ | ETA: 0:01:41
Sampling (1 thread) 42%|█████████████▏ | ETA: 0:01:39
Sampling (1 thread) 43%|█████████████▍ | ETA: 0:01:38
Sampling (1 thread) 44%|█████████████▌ | ETA: 0:01:37
Sampling (1 thread) 44%|█████████████▋ | ETA: 0:01:35
Sampling (1 thread) 44%|█████████████▊ | ETA: 0:01:34
Sampling (1 thread) 45%|██████████████ | ETA: 0:01:33
Sampling (1 thread) 46%|██████████████▏ | ETA: 0:01:31
Sampling (1 thread) 46%|██████████████▎ | ETA: 0:01:30
Sampling (1 thread) 46%|██████████████▍ | ETA: 0:01:29
Sampling (1 thread) 47%|██████████████▋ | ETA: 0:01:28
Sampling (1 thread) 48%|██████████████▊ | ETA: 0:01:26
Sampling (1 thread) 48%|██████████████▉ | ETA: 0:01:25
Sampling (1 thread) 48%|███████████████ | ETA: 0:01:24
Sampling (1 thread) 49%|███████████████▎ | ETA: 0:01:23
Sampling (1 thread) 50%|███████████████▍ | ETA: 0:01:22
Sampling (1 thread) 50%|███████████████▌ | ETA: 0:01:20
┌ Info: Found initial step size
└ ϵ = 0.4
Sampling (1 thread) 50%|███████████████▋ | ETA: 0:01:19
Sampling (1 thread) 51%|███████████████▊ | ETA: 0:01:18
Sampling (1 thread) 52%|████████████████ | ETA: 0:01:17
Sampling (1 thread) 52%|████████████████▏ | ETA: 0:01:16
Sampling (1 thread) 52%|████████████████▎ | ETA: 0:01:15
Sampling (1 thread) 53%|████████████████▍ | ETA: 0:01:14
Sampling (1 thread) 54%|████████████████▋ | ETA: 0:01:13
Sampling (1 thread) 54%|████████████████▊ | ETA: 0:01:11
Sampling (1 thread) 55%|████████████████▉ | ETA: 0:01:10
Sampling (1 thread) 55%|█████████████████ | ETA: 0:01:09
Sampling (1 thread) 56%|█████████████████▎ | ETA: 0:01:08
Sampling (1 thread) 56%|█████████████████▍ | ETA: 0:01:07
Sampling (1 thread) 56%|█████████████████▌ | ETA: 0:01:06
Sampling (1 thread) 57%|█████████████████▋ | ETA: 0:01:05
Sampling (1 thread) 57%|█████████████████▉ | ETA: 0:01:04
Sampling (1 thread) 58%|██████████████████ | ETA: 0:01:03
Sampling (1 thread) 58%|██████████████████▏ | ETA: 0:01:02
Sampling (1 thread) 59%|██████████████████▎ | ETA: 0:01:01
Sampling (1 thread) 60%|██████████████████▌ | ETA: 0:01:00
Sampling (1 thread) 60%|██████████████████▋ | ETA: 0:00:59
Sampling (1 thread) 60%|██████████████████▊ | ETA: 0:00:58
Sampling (1 thread) 61%|██████████████████▉ | ETA: 0:00:58
Sampling (1 thread) 62%|███████████████████▏ | ETA: 0:00:57
Sampling (1 thread) 62%|███████████████████▎ | ETA: 0:00:56
Sampling (1 thread) 62%|███████████████████▍ | ETA: 0:00:55
Sampling (1 thread) 63%|███████████████████▌ | ETA: 0:00:54
Sampling (1 thread) 64%|███████████████████▋ | ETA: 0:00:53
Sampling (1 thread) 64%|███████████████████▉ | ETA: 0:00:52
Sampling (1 thread) 64%|████████████████████ | ETA: 0:00:51
Sampling (1 thread) 65%|████████████████████▏ | ETA: 0:00:50
Sampling (1 thread) 66%|████████████████████▎ | ETA: 0:00:50
Sampling (1 thread) 66%|████████████████████▌ | ETA: 0:00:49
Sampling (1 thread) 66%|████████████████████▋ | ETA: 0:00:48
Sampling (1 thread) 67%|████████████████████▊ | ETA: 0:00:47
Sampling (1 thread) 68%|████████████████████▉ | ETA: 0:00:46
Sampling (1 thread) 68%|█████████████████████▏ | ETA: 0:00:45
Sampling (1 thread) 68%|█████████████████████▎ | ETA: 0:00:45
Sampling (1 thread) 69%|█████████████████████▍ | ETA: 0:00:44
Sampling (1 thread) 70%|█████████████████████▌ | ETA: 0:00:43
Sampling (1 thread) 70%|█████████████████████▊ | ETA: 0:00:42
Sampling (1 thread) 70%|█████████████████████▉ | ETA: 0:00:41
Sampling (1 thread) 71%|██████████████████████ | ETA: 0:00:41
Sampling (1 thread) 72%|██████████████████████▏ | ETA: 0:00:40
Sampling (1 thread) 72%|██████████████████████▍ | ETA: 0:00:39
Sampling (1 thread) 72%|██████████████████████▌ | ETA: 0:00:38
Sampling (1 thread) 73%|██████████████████████▋ | ETA: 0:00:37
Sampling (1 thread) 74%|██████████████████████▊ | ETA: 0:00:37
Sampling (1 thread) 74%|███████████████████████ | ETA: 0:00:36
Sampling (1 thread) 74%|███████████████████████▏ | ETA: 0:00:35
Sampling (1 thread) 75%|███████████████████████▎ | ETA: 0:00:34
┌ Info: Found initial step size
└ ϵ = 0.4
Sampling (1 thread) 76%|███████████████████████▍ | ETA: 0:00:34
Sampling (1 thread) 76%|███████████████████████▌ | ETA: 0:00:33
Sampling (1 thread) 76%|███████████████████████▊ | ETA: 0:00:32
Sampling (1 thread) 77%|███████████████████████▉ | ETA: 0:00:31
Sampling (1 thread) 78%|████████████████████████ | ETA: 0:00:31
Sampling (1 thread) 78%|████████████████████████▏ | ETA: 0:00:30
Sampling (1 thread) 78%|████████████████████████▍ | ETA: 0:00:29
Sampling (1 thread) 79%|████████████████████████▌ | ETA: 0:00:28
Sampling (1 thread) 80%|████████████████████████▋ | ETA: 0:00:28
Sampling (1 thread) 80%|████████████████████████▊ | ETA: 0:00:27
Sampling (1 thread) 80%|█████████████████████████ | ETA: 0:00:26
Sampling (1 thread) 81%|█████████████████████████▏ | ETA: 0:00:25
Sampling (1 thread) 82%|█████████████████████████▎ | ETA: 0:00:25
Sampling (1 thread) 82%|█████████████████████████▍ | ETA: 0:00:24
Sampling (1 thread) 82%|█████████████████████████▋ | ETA: 0:00:23
Sampling (1 thread) 83%|█████████████████████████▊ | ETA: 0:00:23
Sampling (1 thread) 84%|█████████████████████████▉ | ETA: 0:00:22
Sampling (1 thread) 84%|██████████████████████████ | ETA: 0:00:21
Sampling (1 thread) 84%|██████████████████████████▎ | ETA: 0:00:20
Sampling (1 thread) 85%|██████████████████████████▍ | ETA: 0:00:20
Sampling (1 thread) 86%|██████████████████████████▌ | ETA: 0:00:19
Sampling (1 thread) 86%|██████████████████████████▋ | ETA: 0:00:18
Sampling (1 thread) 86%|██████████████████████████▉ | ETA: 0:00:18
Sampling (1 thread) 87%|███████████████████████████ | ETA: 0:00:17
Sampling (1 thread) 88%|███████████████████████████▏ | ETA: 0:00:16
Sampling (1 thread) 88%|███████████████████████████▎ | ETA: 0:00:16
Sampling (1 thread) 88%|███████████████████████████▍ | ETA: 0:00:15
Sampling (1 thread) 89%|███████████████████████████▋ | ETA: 0:00:14
Sampling (1 thread) 90%|███████████████████████████▊ | ETA: 0:00:14
Sampling (1 thread) 90%|███████████████████████████▉ | ETA: 0:00:13
Sampling (1 thread) 90%|████████████████████████████ | ETA: 0:00:12
Sampling (1 thread) 91%|████████████████████████████▎ | ETA: 0:00:12
Sampling (1 thread) 92%|████████████████████████████▍ | ETA: 0:00:11
Sampling (1 thread) 92%|████████████████████████████▌ | ETA: 0:00:10
Sampling (1 thread) 92%|████████████████████████████▋ | ETA: 0:00:10
Sampling (1 thread) 93%|████████████████████████████▉ | ETA: 0:00:09
Sampling (1 thread) 94%|█████████████████████████████ | ETA: 0:00:08
Sampling (1 thread) 94%|█████████████████████████████▏ | ETA: 0:00:08
Sampling (1 thread) 94%|█████████████████████████████▎ | ETA: 0:00:07
Sampling (1 thread) 95%|█████████████████████████████▌ | ETA: 0:00:06
Sampling (1 thread) 96%|█████████████████████████████▋ | ETA: 0:00:06
Sampling (1 thread) 96%|█████████████████████████████▊ | ETA: 0:00:05
Sampling (1 thread) 96%|█████████████████████████████▉ | ETA: 0:00:04
Sampling (1 thread) 97%|██████████████████████████████▏| ETA: 0:00:04
Sampling (1 thread) 98%|██████████████████████████████▎| ETA: 0:00:03
Sampling (1 thread) 98%|██████████████████████████████▍| ETA: 0:00:03
Sampling (1 thread) 98%|██████████████████████████████▌| ETA: 0:00:02
Sampling (1 thread) 99%|██████████████████████████████▊| ETA: 0:00:01
Sampling (1 thread) 100%|██████████████████████████████▉| ETA: 0:00:01
Sampling (1 thread) 100%|███████████████████████████████| Time: 0:02:05
Sampling (1 thread) 100%|███████████████████████████████| Time: 0:02:06
Effective Sample Size is low for venv_N: 398.4840662857016
Effective Sample Size is low for vsys: 281.3710377559386
Rhat is high for vsys: 1.010497541219949
Effective Sample Size is low for a_f: 379.08662726698816
Rhat is high for a_f: 1.0102147830055774
Effective Sample Size is low for m2_i: 272.91431455268133
Rhat is high for m2_i: 1.014274815361154
Effective Sample Size is low for m2_f: 267.0776479641809
Rhat is high for m2_f: 1.0138369601620822
Effective Sample Size is low for vsys_r: 279.84924681059283
Rhat is high for vsys_r: 1.0131153288750523
Rhat is high for K2: 1.0102030886771967
Effective Sample Size is low for vsys_N: 347.0703986637869
Rhat is high for vsys_E: 1.0122591882777747
Rhat is high for m1_f: 1.0139704677873664
Rhat is high for ω_i: 1.010469496463354
Effective Sample Size is low for a_i: 364.50841201186523
Rhat is high for a_i: 1.0129102300144681
Effective Sample Size is low for P_i: 371.8411504649227
Rhat is high for P_i: 1.0138852764089363
Effective Sample Size is low for dm2: 350.52858846535344
Rhat is high for dm2: 1.0111078557158573
Effective Sample Size is low for i_f: 343.3578819091772
Rhat is high for i_f: 1.0100201600733583
Effective Sample Size is low for frac: 280.01413773746435
Rhat is high for frac: 1.0141086923894884
Rhat is high for m1_i: 1.0139704677873664
Effective Sample Size is low for venv_r: 262.39845688513344
Rhat is high for venv_r: 1.0122247853239574
Effective Sample Size is low for vkick: 371.1263186625632Results 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.