Priors
All parameters to your model must have a prior defined. You may provide any continuous, univariate distribution from the Distributions.jl. A few useful distributions include:
NormalUniformLogNormalLogUniformTrucatedNormalVonMises
This pacakge also defines the Sine() distribution for e.g. inclination priors and UniformCircular() for periodic variables. Internally, UniformCircular() creates two standard normal variables and finds the angle between them using arctan. This allows the sampler to smoothly cycle past the ends of the domain. You can specify a different circular domain than (0,2pi) by passing the size of the domain e.g. τ = UniformCircular(1.0).
The VonMise distribution is notable but not commonly used. It is the analog of a normal distribution defined on a circular domain (-π, +π). If you have a Gaussian prior on an angular parameter, a Von Mises distribution is probably more appropriate.
Behind the scenes, Octofitter remaps your parameters to unconstrained domains using the Bijectors.jl (and corrects the priors accordingly). This is essential for good sampling efficiency with HMC based samplers.
This means that e.g. if you define the eccentricity prior as e=Uniform(0,0.5), the sampler will actually generate values across the whole real line and transform them back into the [0,0.5] range before evaluating the orbit. It is therefore essential that your priors do not include invalid domains.
For example, setting a=Normal(3,2) will result in poor sampling efficiency as sometimes negative values for semi-major axis will be drawn (especially if you're using the parallel tempered sampler).
Instead, for parameters like semi-major axis, eccentricity, parallax, and masses, you should truncate any distributions that have negative tails. This can easily be accomplished with TrauncatedNormal or Trunacted(dist, low, high) for any arbitrary distribution.
Kernel Density Estimate Priors
Octofitter has support for sampling from smoothed kernel density estimate priors. These are non-parametric distributions fit to a 1D dataset consisting of random draws. This is one way to include the output of a different model as the prior to a new model. That said, it's usually best to try and incorporate the model directly into the code. There are a few examples on GitHub of this, including atmosphere model grids, cooling tracks, etc.
Using a KDE
First, we will generate some data. In the real world, you would load this data eg. from a CSV file.
using Octofitter, Distributions
# create a smoothed KDE estimate of the samples from a 10+-1 gaussian
kde = Octofitter.KDEDist(randn(1000).+10)KDEDist kernel density estimate distribution
Note that in Octofitter the KDE will have its support truncated to the minimum and maximum values that occur in your dataset, ie. it doesn't allow for infinite long tails.
Now add it to your model as a prior:
planet_b = Planet(
name="b",
basis=Visual{KepOrbit},
likelihoods=[],
variables=@variables begin
a ~ kde # Sample from the KDE here
e ~ Uniform(0.0, 0.99)
i ~ Sine()
M = system.M
ω ~ UniformCircular()
Ω ~ UniformCircular()
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(θ, 50000; M, e, a, i, ω, Ω)
end
)
sys = System(
name="Tutoria",
companions=[planet_b],
likelihoods=[],
variables=@variables begin
M ~ truncated(Normal(1.2, 0.1), lower=0.1)
plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
end
)
model = Octofitter.LogDensityModel(sys)
chain = octofit(model)Chains MCMC chain (1000×29×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.75 seconds
Compute duration = 0.75 seconds
parameters = M, plx, b_a, b_e, b_i, b_ωx, b_ωy, b_Ωx, b_Ωy, b_θx, b_θy, b_ω, b_Ω, b_θ, b_M, b_tp
internals = n_steps, is_accept, acceptance_rate, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, is_adapt, loglike, logpost, tree_depth, numerical_error
Summary Statistics
parameters mean std mcse ess_bulk ess_tail r ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Floa ⋯
M 1.2032 0.1034 0.0036 839.9493 550.5438 0.9 ⋯
plx 49.9986 0.0196 0.0007 845.7689 608.5112 1.0 ⋯
b_a 10.0311 1.0784 0.0329 1089.6567 349.7724 1.0 ⋯
b_e 0.5068 0.2871 0.0099 729.6395 359.5738 1.0 ⋯
b_i 1.5383 0.6715 0.0214 941.3463 702.8012 1.0 ⋯
b_ωx -0.0112 0.7108 0.0352 439.4023 720.8847 1.0 ⋯
b_ωy -0.0196 0.7214 0.0440 288.2921 757.8167 0.9 ⋯
b_Ωx -0.0495 0.7181 0.0427 296.5901 633.5273 1.0 ⋯
b_Ωy -0.0788 0.7021 0.0361 403.4391 850.2053 1.0 ⋯
b_θx -0.0017 0.7004 0.0383 355.1467 687.9573 1.0 ⋯
b_θy 0.0333 0.7329 0.0395 382.7063 738.6797 0.9 ⋯
b_ω 0.0061 1.8210 0.0951 424.4454 679.5607 1.0 ⋯
b_Ω -0.2495 1.8647 0.0826 610.2832 686.4671 1.0 ⋯
b_θ 0.0793 1.7907 0.0848 543.8003 752.2815 1.0 ⋯
b_M 1.2032 0.1034 0.0036 839.9493 550.5438 0.9 ⋯
b_tp 44577.8373 4217.1820 149.8135 859.8357 799.3961 0.9 ⋯
2 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
M 1.0036 1.1360 1.2004 1.2715 1.4099
plx 49.9614 49.9852 49.9992 50.0112 50.0370
b_a 8.0080 9.2891 10.0166 10.7825 12.0118
b_e 0.0169 0.2627 0.5163 0.7584 0.9663
b_i 0.3201 1.0239 1.5349 2.0352 2.7911
b_ωx -1.0807 -0.7118 -0.0036 0.6734 1.0674
b_ωy -1.0602 -0.7308 -0.0144 0.6914 1.0788
b_Ωx -1.0628 -0.7467 -0.1158 0.6630 1.0459
b_Ωy -1.0633 -0.7421 -0.1546 0.5978 1.0510
b_θx -1.0598 -0.6634 -0.0026 0.6677 1.0887
b_θy -1.0585 -0.7012 0.0726 0.7668 1.0611
b_ω -2.9499 -1.5592 -0.0404 1.6231 3.0122
b_Ω -3.0457 -1.8875 -0.4049 1.3877 2.9854
b_θ -2.9633 -1.4932 0.1215 1.6175 2.9381
b_M 1.0036 1.1360 1.2004 1.2715 1.4099
b_tp 37620.4104 40719.2805 44547.2410 48959.9744 49985.6974
We now examine the posterior and verify that it matches our KDE prior:
dat = chain[:b_a][:]
@show mean(dat) std(dat)mean(dat) = 10.031054759955724
std(dat) = 1.0783668738628527Observable Based Priors
Octofitter implements observable-based priors from O'Neil 2019 for relative astrometry. You can fit a model to astrometry using observable-based priors using the following recipe:
using Octofitter, Distributions
astrom_dat = Table(
epoch=[mjd("2020-12-20")],
ra=[400.0],
σ_ra=[5.0],
dec=[400.0],
σ_dec=[5.0]
)
astrom_like = PlanetRelAstromLikelihood(astrom_dat, name="rel astrom. 1")
planet_b = Planet(
name="b",
basis=Visual{KepOrbit},
likelihoods=[ObsPriorAstromONeil2019(astrom_like)],
variables=@variables begin
# For using with ObsPriors:
P ~ Uniform(0.001, 1000)
M = system.M
a = cbrt(M * P^2)
e ~ Uniform(0.0, 1.0)
i ~ Sine()
ω ~ UniformCircular()
Ω ~ UniformCircular()
mass ~ LogUniform(0.01, 100)
τ ~ UniformCircular(1.0)
tp = τ*P*365.25 + 58849 # reference epoch for τ. Choose an MJD date near your data.
end
)
sys = System(
name="System1",
companions=[planet_b],
likelihoods=[],
variables=@variables begin
plx ~ Normal(21.219, 0.060)
M ~ truncated(Normal(1.1, 0.2),lower=0.1)
end
)