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 Distributions

We 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/mcqES/src/sample.jl:432
Sampling (1 thread)   0%|                               |  ETA: N/A
┌ Info: Found initial step size
└   ϵ = 0.4
Sampling (1 thread)   0%|▏                              |  ETA: 3:10:15
Sampling (1 thread)   1%|▎                              |  ETA: 1:35:01
Sampling (1 thread)   2%|▌                              |  ETA: 1:03:28
Sampling (1 thread)   2%|▋                              |  ETA: 0:47:41
Sampling (1 thread)   2%|▊                              |  ETA: 0:38:13
Sampling (1 thread)   3%|▉                              |  ETA: 0:31:54
Sampling (1 thread)   4%|█▏                             |  ETA: 0:27:23
Sampling (1 thread)   4%|█▎                             |  ETA: 0:24:00
Sampling (1 thread)   4%|█▍                             |  ETA: 0:21:21
Sampling (1 thread)   5%|█▌                             |  ETA: 0:19:14
Sampling (1 thread)   6%|█▊                             |  ETA: 0:17:31
Sampling (1 thread)   6%|█▉                             |  ETA: 0:16:04
Sampling (1 thread)   6%|██                             |  ETA: 0:14:50
Sampling (1 thread)   7%|██▏                            |  ETA: 0:13:48
Sampling (1 thread)   8%|██▍                            |  ETA: 0:12:53
Sampling (1 thread)   8%|██▌                            |  ETA: 0:12:05
Sampling (1 thread)   8%|██▋                            |  ETA: 0:11:24
Sampling (1 thread)   9%|██▊                            |  ETA: 0:10:47
Sampling (1 thread)  10%|███                            |  ETA: 0:10:14
Sampling (1 thread)  10%|███▏                           |  ETA: 0:09:43
Sampling (1 thread)  10%|███▎                           |  ETA: 0:09:16
Sampling (1 thread)  11%|███▍                           |  ETA: 0:08:51
Sampling (1 thread)  12%|███▋                           |  ETA: 0:08:28
Sampling (1 thread)  12%|███▊                           |  ETA: 0:08:07
Sampling (1 thread)  12%|███▉                           |  ETA: 0:07:48
Sampling (1 thread)  13%|████                           |  ETA: 0:07:30
Sampling (1 thread)  14%|████▏                          |  ETA: 0:07:13
Sampling (1 thread)  14%|████▍                          |  ETA: 0:06:58
Sampling (1 thread)  14%|████▌                          |  ETA: 0:06:44
Sampling (1 thread)  15%|████▋                          |  ETA: 0:06:30
Sampling (1 thread)  16%|████▊                          |  ETA: 0:06:17
Sampling (1 thread)  16%|█████                          |  ETA: 0:06:06
Sampling (1 thread)  16%|█████▏                         |  ETA: 0:05:54
Sampling (1 thread)  17%|█████▎                         |  ETA: 0:05:44
Sampling (1 thread)  18%|█████▍                         |  ETA: 0:05:34
Sampling (1 thread)  18%|█████▋                         |  ETA: 0:05:25
Sampling (1 thread)  18%|█████▊                         |  ETA: 0:05:16
Sampling (1 thread)  19%|█████▉                         |  ETA: 0:05:07
Sampling (1 thread)  20%|██████                         |  ETA: 0:04:59
Sampling (1 thread)  20%|██████▎                        |  ETA: 0:04:51
Sampling (1 thread)  20%|██████▍                        |  ETA: 0:04:44
Sampling (1 thread)  21%|██████▌                        |  ETA: 0:04:37
Sampling (1 thread)  22%|██████▋                        |  ETA: 0:04:30
Sampling (1 thread)  22%|██████▉                        |  ETA: 0:04:24
Sampling (1 thread)  22%|███████                        |  ETA: 0:04:18
Sampling (1 thread)  23%|███████▏                       |  ETA: 0:04:12
Sampling (1 thread)  24%|███████▎                       |  ETA: 0:04:06
Sampling (1 thread)  24%|███████▌                       |  ETA: 0:04:01
Sampling (1 thread)  24%|███████▋                       |  ETA: 0:03:56
Sampling (1 thread)  25%|███████▊                       |  ETA: 0:03:51
┌ Info: Found initial step size
└   ϵ = 0.4
Sampling (1 thread)  26%|███████▉                       |  ETA: 0:04:14
Sampling (1 thread)  26%|████████                       |  ETA: 0:04:08
Sampling (1 thread)  26%|████████▎                      |  ETA: 0:04:03
Sampling (1 thread)  27%|████████▍                      |  ETA: 0:03:58
Sampling (1 thread)  28%|████████▌                      |  ETA: 0:03:53
Sampling (1 thread)  28%|████████▋                      |  ETA: 0:03:48
Sampling (1 thread)  28%|████████▉                      |  ETA: 0:03:44
Sampling (1 thread)  29%|█████████                      |  ETA: 0:03:39
Sampling (1 thread)  30%|█████████▏                     |  ETA: 0:03:35
Sampling (1 thread)  30%|█████████▎                     |  ETA: 0:03:31
Sampling (1 thread)  30%|█████████▌                     |  ETA: 0:03:27
Sampling (1 thread)  31%|█████████▋                     |  ETA: 0:03:23
Sampling (1 thread)  32%|█████████▊                     |  ETA: 0:03:19
Sampling (1 thread)  32%|█████████▉                     |  ETA: 0:03:15
Sampling (1 thread)  32%|██████████▏                    |  ETA: 0:03:11
Sampling (1 thread)  33%|██████████▎                    |  ETA: 0:03:07
Sampling (1 thread)  34%|██████████▍                    |  ETA: 0:03:04
Sampling (1 thread)  34%|██████████▌                    |  ETA: 0:03:01
Sampling (1 thread)  34%|██████████▊                    |  ETA: 0:02:57
Sampling (1 thread)  35%|██████████▉                    |  ETA: 0:02:54
Sampling (1 thread)  36%|███████████                    |  ETA: 0:02:51
Sampling (1 thread)  36%|███████████▏                   |  ETA: 0:02:48
Sampling (1 thread)  36%|███████████▍                   |  ETA: 0:02:45
Sampling (1 thread)  37%|███████████▌                   |  ETA: 0:02:42
Sampling (1 thread)  38%|███████████▋                   |  ETA: 0:02:40
Sampling (1 thread)  38%|███████████▊                   |  ETA: 0:02:37
Sampling (1 thread)  38%|███████████▉                   |  ETA: 0:02:34
Sampling (1 thread)  39%|████████████▏                  |  ETA: 0:02:32
Sampling (1 thread)  40%|████████████▎                  |  ETA: 0:02:29
Sampling (1 thread)  40%|████████████▍                  |  ETA: 0:02:27
Sampling (1 thread)  40%|████████████▌                  |  ETA: 0:02:24
Sampling (1 thread)  41%|████████████▊                  |  ETA: 0:02:22
Sampling (1 thread)  42%|████████████▉                  |  ETA: 0:02:19
Sampling (1 thread)  42%|█████████████                  |  ETA: 0:02:17
Sampling (1 thread)  42%|█████████████▏                 |  ETA: 0:02:15
Sampling (1 thread)  43%|█████████████▍                 |  ETA: 0:02:13
Sampling (1 thread)  44%|█████████████▌                 |  ETA: 0:02:11
Sampling (1 thread)  44%|█████████████▋                 |  ETA: 0:02:08
Sampling (1 thread)  44%|█████████████▊                 |  ETA: 0:02:06
Sampling (1 thread)  45%|██████████████                 |  ETA: 0:02:04
Sampling (1 thread)  46%|██████████████▏                |  ETA: 0:02:02
Sampling (1 thread)  46%|██████████████▎                |  ETA: 0:02:00
Sampling (1 thread)  46%|██████████████▍                |  ETA: 0:01:58
Sampling (1 thread)  47%|██████████████▋                |  ETA: 0:01:57
Sampling (1 thread)  48%|██████████████▊                |  ETA: 0:01:55
Sampling (1 thread)  48%|██████████████▉                |  ETA: 0:01:53
Sampling (1 thread)  48%|███████████████                |  ETA: 0:01:51
Sampling (1 thread)  49%|███████████████▎               |  ETA: 0:01:49
Sampling (1 thread)  50%|███████████████▍               |  ETA: 0:01:47
Sampling (1 thread)  50%|███████████████▌               |  ETA: 0:01:46
┌ Info: Found initial step size
└   ϵ = 0.4
Sampling (1 thread)  50%|███████████████▋               |  ETA: 0:01:44
Sampling (1 thread)  51%|███████████████▊               |  ETA: 0:01:42
Sampling (1 thread)  52%|████████████████               |  ETA: 0:01:41
Sampling (1 thread)  52%|████████████████▏              |  ETA: 0:01:39
Sampling (1 thread)  52%|████████████████▎              |  ETA: 0:01:38
Sampling (1 thread)  53%|████████████████▍              |  ETA: 0:01:36
Sampling (1 thread)  54%|████████████████▋              |  ETA: 0:01:34
Sampling (1 thread)  54%|████████████████▊              |  ETA: 0:01:33
Sampling (1 thread)  55%|████████████████▉              |  ETA: 0:01:31
Sampling (1 thread)  55%|█████████████████              |  ETA: 0:01:30
Sampling (1 thread)  56%|█████████████████▎             |  ETA: 0:01:28
Sampling (1 thread)  56%|█████████████████▍             |  ETA: 0:01:27
Sampling (1 thread)  56%|█████████████████▌             |  ETA: 0:01:25
Sampling (1 thread)  57%|█████████████████▋             |  ETA: 0:01:24
Sampling (1 thread)  57%|█████████████████▉             |  ETA: 0:01:23
Sampling (1 thread)  58%|██████████████████             |  ETA: 0:01:21
Sampling (1 thread)  58%|██████████████████▏            |  ETA: 0:01:20
Sampling (1 thread)  59%|██████████████████▎            |  ETA: 0:01:19
Sampling (1 thread)  60%|██████████████████▌            |  ETA: 0:01:17
Sampling (1 thread)  60%|██████████████████▋            |  ETA: 0:01:16
Sampling (1 thread)  60%|██████████████████▊            |  ETA: 0:01:15
Sampling (1 thread)  61%|██████████████████▉            |  ETA: 0:01:13
Sampling (1 thread)  62%|███████████████████▏           |  ETA: 0:01:12
Sampling (1 thread)  62%|███████████████████▎           |  ETA: 0:01:11
Sampling (1 thread)  62%|███████████████████▍           |  ETA: 0:01:09
Sampling (1 thread)  63%|███████████████████▌           |  ETA: 0:01:08
Sampling (1 thread)  64%|███████████████████▋           |  ETA: 0:01:07
Sampling (1 thread)  64%|███████████████████▉           |  ETA: 0:01:06
Sampling (1 thread)  64%|████████████████████           |  ETA: 0:01:05
Sampling (1 thread)  65%|████████████████████▏          |  ETA: 0:01:03
Sampling (1 thread)  66%|████████████████████▎          |  ETA: 0:01:02
Sampling (1 thread)  66%|████████████████████▌          |  ETA: 0:01:01
Sampling (1 thread)  66%|████████████████████▋          |  ETA: 0:01:00
Sampling (1 thread)  67%|████████████████████▊          |  ETA: 0:00:59
Sampling (1 thread)  68%|████████████████████▉          |  ETA: 0:00:58
Sampling (1 thread)  68%|█████████████████████▏         |  ETA: 0:00:57
Sampling (1 thread)  68%|█████████████████████▎         |  ETA: 0:00:56
Sampling (1 thread)  69%|█████████████████████▍         |  ETA: 0:00:54
Sampling (1 thread)  70%|█████████████████████▌         |  ETA: 0:00:53
Sampling (1 thread)  70%|█████████████████████▊         |  ETA: 0:00:52
Sampling (1 thread)  70%|█████████████████████▉         |  ETA: 0:00:51
Sampling (1 thread)  71%|██████████████████████         |  ETA: 0:00:50
Sampling (1 thread)  72%|██████████████████████▏        |  ETA: 0:00:49
Sampling (1 thread)  72%|██████████████████████▍        |  ETA: 0:00:48
Sampling (1 thread)  72%|██████████████████████▌        |  ETA: 0:00:47
Sampling (1 thread)  73%|██████████████████████▋        |  ETA: 0:00:46
Sampling (1 thread)  74%|██████████████████████▊        |  ETA: 0:00:45
Sampling (1 thread)  74%|███████████████████████        |  ETA: 0:00:44
Sampling (1 thread)  74%|███████████████████████▏       |  ETA: 0:00:43
Sampling (1 thread)  75%|███████████████████████▎       |  ETA: 0:00:42
┌ Info: Found initial step size
└   ϵ = 0.4
Sampling (1 thread)  76%|███████████████████████▍       |  ETA: 0:00:41
Sampling (1 thread)  76%|███████████████████████▌       |  ETA: 0:00:40
Sampling (1 thread)  76%|███████████████████████▊       |  ETA: 0:00:39
Sampling (1 thread)  77%|███████████████████████▉       |  ETA: 0:00:38
Sampling (1 thread)  78%|████████████████████████       |  ETA: 0:00:37
Sampling (1 thread)  78%|████████████████████████▏      |  ETA: 0:00:36
Sampling (1 thread)  78%|████████████████████████▍      |  ETA: 0:00:35
Sampling (1 thread)  79%|████████████████████████▌      |  ETA: 0:00:34
Sampling (1 thread)  80%|████████████████████████▋      |  ETA: 0:00:33
Sampling (1 thread)  80%|████████████████████████▊      |  ETA: 0:00:32
Sampling (1 thread)  80%|█████████████████████████      |  ETA: 0:00:32
Sampling (1 thread)  81%|█████████████████████████▏     |  ETA: 0:00:31
Sampling (1 thread)  82%|█████████████████████████▎     |  ETA: 0:00:30
Sampling (1 thread)  82%|█████████████████████████▍     |  ETA: 0:00:29
Sampling (1 thread)  82%|█████████████████████████▋     |  ETA: 0:00:28
Sampling (1 thread)  83%|█████████████████████████▊     |  ETA: 0:00:27
Sampling (1 thread)  84%|█████████████████████████▉     |  ETA: 0:00:26
Sampling (1 thread)  84%|██████████████████████████     |  ETA: 0:00:25
Sampling (1 thread)  84%|██████████████████████████▎    |  ETA: 0:00:24
Sampling (1 thread)  85%|██████████████████████████▍    |  ETA: 0:00:24
Sampling (1 thread)  86%|██████████████████████████▌    |  ETA: 0:00:23
Sampling (1 thread)  86%|██████████████████████████▋    |  ETA: 0:00:22
Sampling (1 thread)  86%|██████████████████████████▉    |  ETA: 0:00:21
Sampling (1 thread)  87%|███████████████████████████    |  ETA: 0:00:20
Sampling (1 thread)  88%|███████████████████████████▏   |  ETA: 0:00:19
Sampling (1 thread)  88%|███████████████████████████▎   |  ETA: 0:00:19
Sampling (1 thread)  88%|███████████████████████████▍   |  ETA: 0:00:18
Sampling (1 thread)  89%|███████████████████████████▋   |  ETA: 0:00:17
Sampling (1 thread)  90%|███████████████████████████▊   |  ETA: 0:00:16
Sampling (1 thread)  90%|███████████████████████████▉   |  ETA: 0:00:15
Sampling (1 thread)  90%|████████████████████████████   |  ETA: 0:00:14
Sampling (1 thread)  91%|████████████████████████████▎  |  ETA: 0:00:14
Sampling (1 thread)  92%|████████████████████████████▍  |  ETA: 0:00:13
Sampling (1 thread)  92%|████████████████████████████▌  |  ETA: 0:00:12
Sampling (1 thread)  92%|████████████████████████████▋  |  ETA: 0:00:11
Sampling (1 thread)  93%|████████████████████████████▉  |  ETA: 0:00:11
Sampling (1 thread)  94%|█████████████████████████████  |  ETA: 0:00:10
Sampling (1 thread)  94%|█████████████████████████████▏ |  ETA: 0:00:09
Sampling (1 thread)  94%|█████████████████████████████▎ |  ETA: 0:00:08
Sampling (1 thread)  95%|█████████████████████████████▌ |  ETA: 0:00:07
Sampling (1 thread)  96%|█████████████████████████████▋ |  ETA: 0:00:07
Sampling (1 thread)  96%|█████████████████████████████▊ |  ETA: 0:00:06
Sampling (1 thread)  96%|█████████████████████████████▉ |  ETA: 0:00:05
Sampling (1 thread)  97%|██████████████████████████████▏|  ETA: 0:00:04
Sampling (1 thread)  98%|██████████████████████████████▎|  ETA: 0:00:04
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:24
Sampling (1 thread) 100%|███████████████████████████████| Time: 0:02:25
Rhat is high for venv_N: 1.0106549286617872
Effective Sample Size is low for vsys: 185.78305230013757
Rhat is high for vsys: 1.0154253065621277
Effective Sample Size is low for m2_i: 306.39540264630307
Effective Sample Size is low for m2_f: 332.4299655767324
Effective Sample Size is low for vsys_r: 249.6455839544729
Rhat is high for vsys_r: 1.0154573997380727
Effective Sample Size is low for K2: 398.7104852513822
Effective Sample Size is low for vsys_E: 382.6163026254078
Rhat is high for m1_f: 1.0118064389783574
Effective Sample Size is low for venv_E: 358.1691829524051
Rhat is high for venv_E: 1.0105030571579225
Effective Sample Size is low for a_i: 380.6066635075972
Rhat is high for a_i: 1.0104139648551551
Effective Sample Size is low for P_i: 214.9697423857309
Rhat is high for P_i: 1.0210214882829436
Effective Sample Size is low for dm2: 273.693870562723
Effective Sample Size is low for frac: 196.73453045963484
Rhat is high for frac: 1.0167647813618035
Rhat is high for m1_i: 1.0118064389783574
Effective Sample Size is low for venv_r: 259.7113713896337
Rhat is high for venv_r: 1.0142737146433387
Effective Sample Size is low for vkick: 222.31938070861528

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.