Multi-Planet RV Fits

This tutorial shows how to perform a multi-planet RV fit, and compare the bayesian evidence between the two models.

using Octofitter
using OctofitterRadialVelocity
using CairoMakie
using PairPlots
using Distributions
using PlanetOrbits

To begin, we create simulated data. We imagine that we have two different instruments.

using Random
Random.seed!(1)

orb_template_1 = orbit(a = 1.0,e = 0.05,ω = 1π/4,M = 1.0,tp =58800)
mass_1 = 0.25*1e-3
orb_template_2 = orbit(a = 5.0,e = 0.4,ω = 1π/4,M = 1.0,tp =59800)
mass_2 = 1.0*1e-3

epochs = (58400:150:69400) .+ 10 .* randn.()
rv = radvel.(orb_template_1, epochs, mass_1) .+ radvel.(orb_template_2, epochs, mass_2)
rvlike1 = MarginalizedStarAbsoluteRVLikelihood(
    Table(epoch=epochs, rv=rv .+ 4 .* randn.(), σ_rv=[4 .* abs.(randn.()) .+ 1 for _ in 1:length(epochs)]),
    jitter=:jitter1,
    instrument_name="DATA 1"
)

epochs = (65400:100:71400) .+ 10 .* randn.()
rv = radvel.(orb_template_1, epochs, mass_1) .+ radvel.(orb_template_2, epochs, mass_2)
rvlike2 = MarginalizedStarAbsoluteRVLikelihood(
    Table(epoch=epochs, rv=rv .+ 2 .* randn.() .+ 7, σ_rv=[2 .* abs.(randn.()) .+ 1 for _ in 1:length(epochs)]),
    jitter=:jitter2,
    instrument_name="DATA 2"
)

fig = Figure()
ax = Axis(
    fig[1,1],
    xlabel="epoch [mjd]",
    ylabel="rv [m/s]"
)
Makie.scatter!(ax, rvlike1.table.epoch, rvlike1.table.rv)
Makie.errorbars!(ax, rvlike1.table.epoch, rvlike1.table.rv, rvlike1.table.σ_rv)
Makie.scatter!(ax, rvlike2.table.epoch, rvlike2.table.rv)
Makie.errorbars!(ax, rvlike2.table.epoch, rvlike2.table.rv, rvlike2.table.σ_rv)
fig

Two Planet Model

@planet b RadialVelocityOrbit begin
    e ~ Uniform(0,0.999999)
    mass ~ Uniform(0, 10)
    i ~ Sine()
    ω ~ Uniform(0,2pi)
    τ ~ Uniform(0,1.0)

    P_kep_yrs ~ Uniform(0, 100)
    a = ∛(system.M * b.P_kep_yrs^2)
    tp =  b.τ*b.P_kep_yrs*365.25 + 58400
end

@planet c RadialVelocityOrbit begin
    e ~ Uniform(0,0.999999)
    mass ~ Uniform(0, 10)
    i ~ Sine()
    ω ~ Uniform(0,2pi)
    τ ~ Uniform(0,1.0)

    P_kep_yrs ~ Uniform(0, 100)
    a = ∛(system.M * c.P_kep_yrs^2)
    tp =  c.τ*c.P_kep_yrs*365.25 + 58400
end

@system sim_2p begin
    M = 1.0
    jitter1 ~ Uniform(0, 20_000)
    jitter2 ~ Uniform(0, 20_000)
end rvlike1 rvlike2 b c

model_2p = Octofitter.LogDensityModel(sim_2p)

LogDensityModel for System sim_2p of dimension 14 and 135 epochs with fields .ℓπcallback and .∇ℓπcallback

Sample from the posterior

using Pigeons
results_2p, pt_2p = octofit_pigeons(model_2p, n_rounds=10)

(chain = MCMC chain (1024×23×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Plot RV curve, phase folded curve, and binned residuals:

Octofitter.rvpostplot(model_2p, results_2p)

One Planet Model

We now create a new system object that only includes one planet (we dropped c, in this case).

@system sim_1p begin
    M = 1.0 #~ truncated(Normal(1, 0.04),lower=0.1) # (Baines & Armstrong 2011).
    jitter1 ~ Uniform(0, 20_000)
    jitter2 ~ Uniform(0, 20_000)
end rvlike1 rvlike2 b

model_1p = Octofitter.LogDensityModel(sim_1p)

LogDensityModel for System sim_1p of dimension 8 and 135 epochs with fields .ℓπcallback and .∇ℓπcallback

Sample from the posterior

using Pigeons
results_1p, pt_1p = octofit_pigeons(model_1p, n_rounds=10)

(chain = MCMC chain (1024×15×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

Plot RV curve, phase folded curve, and binned residuals:

Octofitter.rvpostplot(model_1p, results_1p)

Model Comparison: Bayesian Evidence

Octofitter with Pigeons directly calculates the log Bayesian evidence using the "stepping stone" method. This should be more reliable than even nested sampling, and certainly more reliable than approximate methods like the BIC etc.

Z1 = stepping_stone(pt_1p)
Z2 = stepping_stone(pt_2p)

log_BF₁₀ = Z2-Z1

121.48007246948146

Here is a standard guideline you can use to interpret the evidence:

As you can see, the evidence for there being two planets is "extreme" in this case. Try adjusting the masses of the two planets and see how this changes!

Parameterizations

When using the evidence for model comparisons, a model with more specific priors will have more evidence than an quivalent model with broad priors.

In our two planet model above, we made two exactly equivalent planets. If you inspect the chains, you may notice that the two planets often flip back and forth – sometimes b has the longer period, and sometimes c does.

For example, here is a histogram of the period of planet b:

hist(vec(results_2p[:b_P_kep_yrs]), bins=100)

We can refine the two planet model a bit by adjusting the priors such that planet c always has a longer period than planet b.

This will make analysis a little more straightforward, but crucially it will also increase the evidence of this model, by approximately halving the prior volume–-thus making a more specific prediction.

There are several ways we could do this. Here, we add a "nominal period" variable and reparameterize the two planets as ratios of this nominal period.

@planet b RadialVelocityOrbit begin
    e ~ Uniform(0,0.999999)
    mass ~ Uniform(0, 10)
    i ~ Sine()
    ω ~ Uniform(0,2pi)
    τ ~ Uniform(0,1.0)

    P_kep_yrs = system.P_yrs_nom * system.P_ratio_b
    a = ∛(system.M * b.P_kep_yrs^2)
    tp =  b.τ*b.P_kep_yrs*365.25 + 58400
end

@planet c RadialVelocityOrbit begin
    e ~ Uniform(0,0.999999)
    mass ~ Uniform(0, 10)
    i ~ Sine()
    ω ~ Uniform(0,2pi)
    τ ~ Uniform(0,1.0)

    P_kep_yrs = system.P_yrs_nom * system.P_ratio_c
    a = ∛(system.M * c.P_kep_yrs^2)
    tp =  c.τ*c.P_kep_yrs*365.25 + 58400
end

@system sim_2p_v2 begin
    M = 1.0
    jitter1 ~ Uniform(0, 20_000)
    jitter2 ~ Uniform(0, 20_000)

    P_yrs_nom ~ Uniform(0, 100)
    P_ratio_b ~ Uniform(0, 0.5)
    P_ratio_c ~ Uniform(0.5, 1)
end rvlike1 rvlike2 b c

model_2p_v2 = Octofitter.LogDensityModel(sim_2p_v2)

LogDensityModel for System sim_2p_v2 of dimension 15 and 135 epochs with fields .ℓπcallback and .∇ℓπcallback

Sample from the posterior

using Pigeons
results_2p_v2, pt_2p_v2 = octofit_pigeons(model_2p_v2, n_rounds=10)

(chain = MCMC chain (1024×26×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))

The planet with the wider orbit is now consistently plotted in the bottom panel (meaning that planet b and c are no longer trading back and forth):

Octofitter.rvpostplot(model_2p_v2, results_2p_v2)

If we look again at the log-evidence, we see that this parameterization (Z3) is even more favoured. This is because this small change in parameterization makes considerably

Z1 = stepping_stone(pt_1p)
Z2 = stepping_stone(pt_2p)
Z3 = stepping_stone(pt_2p_v2)

Z1, Z2, Z3

(-917.2230277100025, -795.742955240521, -793.607352029485)

As a final treat, let's animate the orbit plots. All the previous images were visualizing a single posterior draw. In this animation, we'll loop over many different samples:

Octofitter.rvpostplot_animated(model_2p_v2, results_2p_v2)

"rv-posterior.mp4"

prior_model = Octofitter.LogDensityModel(Octofitter.prior_only_model(model_1p.system, exclude_all=true))
_, pt_prior = octofit_pigeons(prior_model, n_rounds=10) # should be very quick!
log_Z0 = stepping_stone(pt_prior)

-0.004672803525865188

log_Z1_over_Z0 = stepping_stone(pt_1p)
log_Z1 = log_Z1_over_Z0 - log_Z0

-917.2183549064766