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:

  • Normal
  • Uniform
  • LogNormal
  • LogUniform
  • TrucatedNormal
  • VonMises

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 Visual{KepOrbit} begin
    a ~ kde # Sample from the KDE here
    e ~ Uniform(0.0, 0.99)
    i ~ Sine()
    ω ~ UniformCircular()
    Ω ~ UniformCircular()
    θ ~ UniformCircular()
    tp = θ_at_epoch_to_tperi(system,b,50000)
end
@system Tutoria begin
    M ~ truncated(Normal(1.2, 0.1), lower=0.1)
    plx ~ truncated(Normal(50.0, 0.02), lower=0.1)
end b
model = Octofitter.LogDensityModel(Tutoria)
chain = octofit(model)
Chains MCMC chain (1000×28×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 0.73 seconds
Compute duration  = 0.73 seconds
parameters        = M, plx, b_a, b_e, b_i, b_ωy, b_ωx, b_Ωy, b_Ωx, b_θy, b_θx, b_ω, b_Ω, b_θ, 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_ωy      -0.0112      0.7108     0.0352    439.4023   720.8847    1.0 ⋯
        b_ωx      -0.0196      0.7214     0.0440    288.2921   757.8167    0.9 ⋯
        b_Ωy      -0.0495      0.7181     0.0427    296.5901   633.5273    1.0 ⋯
        b_Ωx      -0.0788      0.7021     0.0361    403.4391   850.2053    1.0 ⋯
        b_θy      -0.0017      0.7004     0.0383    355.1467   687.9573    1.0 ⋯
        b_θx       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_tp   44997.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_ωy      -1.0807      -0.7118      -0.0036       0.6734       1.0674
        b_ωx      -1.0602      -0.7308      -0.0144       0.6914       1.0788
        b_Ωy      -1.0628      -0.7467      -0.1158       0.6630       1.0459
        b_Ωx      -1.0633      -0.7421      -0.1546       0.5978       1.0510
        b_θy      -1.0598      -0.6634      -0.0026       0.6677       1.0887
        b_θx      -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_tp   38040.4104   41139.2805   44967.2410   49379.9744   50405.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.0783668738628527

Observable 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_like = PlanetRelAstromLikelihood(
    (;epoch=mjd("2020-12-20"), ra=400.0, σ_ra=5.0, dec=400.0, σ_dec=5.0)
)

@planet b Visual{KepOrbit} begin
    # For using with ObsPriors:
    P ~ Uniform(0.001, 1000)
    a = cbrt(system.M * b.P^2)

    e ~ Uniform(0.0, 1.0)
    i ~ Sine()
    ω ~ UniformCircular()
    Ω ~ UniformCircular()
    mass ~ LogUniform(0.01, 100)

    τ ~ UniformCircular(1.0)
    tp =  b.τ*b.P*365.25 + 58849 # reference epoch for τ. Choose an MJD date near your data.
end astrom_like ObsPriorAstromONeil2019(astrom_like);

@system System1 begin
    plx ~ Normal(21.219, 0.060)
	M ~ truncated(Normal(1.1, 0.2),lower=0.1)
end b