The Exoplanet Simple Orbit Fitting Toolbox (ExoSOFT): An Open-source Tool for Efficient Fitting of Astrometric and Radial Velocity Data

Model fitting can be broken down into two main tasks: solving for a set of parameters approaching those with the maximum likelihood and estimating their uncertainties. In a Bayesian framework, the full set of information is contained within the posterior probability distribution of the model parameters, $p(\mathrm{Model}| \mathrm{Data})$. Ford (2005, 2006) discuss the use of Bayesian inference to solve for this distribution in the case of fitting observations of exoplanets to Keplerian models; for details beyond those given here, we refer the reader to those papers.

The posterior probability distribution, $p(\mathrm{Model}| \mathrm{Data})$, is computed using Bayes' theorem,

$$\begin{eqnarray}&&p(\mathrm{Model}| \mathrm{Data})\propto p(\mathrm{Data}| \mathrm{Model})p(\mathrm{Model}).\end{eqnarray} \tag{ 16 }$$

Here $p(\mathrm{Data}| \mathrm{Model})$ is the likelihood of the model parameters, ${\mathscr{L}}(\mathrm{Model})$, and p(Model) is the prior probability of the parameters based on previous knowledge. Assuming that the errors in the observed data are independent and Gaussian distributed, the likelihood can be expressed in terms of the usual χ²:

$$\begin{eqnarray}&&{\mathscr{L}}(\mathrm{Model})=p(\mathrm{Data}| \mathrm{Model})\propto \exp \left(-\displaystyle \frac{{\chi }^{2}}{2}\right),\end{eqnarray} \tag{ 17 }$$

with χ² being a sum over data points i,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\,\displaystyle \frac{{({\mathrm{Model}}_{i}-{\mathrm{Data}}_{i})}^{2}}{{\sigma }_{i}^{2}}\end{eqnarray} \tag{ 18 }$$

The posterior probability distribution can then be computed by finding the value of ${\mathscr{L}}(\mathrm{Model})p(\mathrm{Model})$ for all possible parameter combinations and integrating. Due to the high dimensionality of Keplerian models however, the parameter space is large, making this brute-force approach computationally problematic. A variety of routines exist for sufficiently sampling the posterior without calculating its value everywhere. In Section 3.2 we will discuss the approaches used in ExoSOFT, particularly that of Markov Chain Monte Carlo (MCMC), along with the specific priors for each of the model parameters.

Forming the posterior distributions can be achieved by integrating Bayes' theorem over the model's parameter space. Depending on the form of the likelihood function and the dimensionality of the model space, evaluating this integral directly or through "shotgun" Monte Carlo can prove computationally impractical. One alternative is to draw statistically dependent samples with MCMC (Liu 2004). A Markov Chain is a sequence of points satisfying the dependency defined by the Markov property: each point in the chain is a function only of its immediate predecessor with no memory of the chain history. MCMC guarantees convergence to the posterior probability distributions in the limit of a large number of steps. A widely used algorithm to form such chains is that of Metropolis–Hastings, originally detailed in Metropolis et al. (1953) and further generalized by Hastings (1970). Unfortunately, the convergence rate is typically impossible to compute; we discuss convergence in more detail in Section 3.3.

There are multiple approaches to sampling the posterior distributions of the Keplerian model in Section 2 with ExoSOFT. One is to use a built-in MCMC sampler based on the Metropolis–Hastings algorithm. This implements the standard form of the rejection function,

$$\begin{eqnarray}&&r({\boldsymbol{\xi }},{{\boldsymbol{\xi }}}_{j+1})=\min \,\left\{1,\displaystyle \frac{{\mathscr{L}}({{\boldsymbol{\xi }}}_{j+1})}{{\mathscr{L}}({{\boldsymbol{\xi }}}_{j})}\displaystyle \frac{p({{\boldsymbol{\xi }}}_{j+1})}{p({{\boldsymbol{\xi }}}_{j})}\right\},\end{eqnarray} \tag{ 19 }$$

where ${{\boldsymbol{\xi }}}_{j+1}$ is a step, a set of parameters; with this, the likelihood is ${\mathscr{L}}={\mathscr{L}}({\mathrm{Model}}_{j})={\mathscr{L}}({{\boldsymbol{\xi }}}_{j})$, and the prior probability is $p=p(\mathrm{Model})=p({{\boldsymbol{\xi }}}_{j})$. The proposed steps are drawn from normal distributions of fixed width in each parameter centered on the current step. Following the Gaussian error assumption in Section 3.1, the likelihood ratio is

$$\begin{eqnarray}&&\displaystyle \frac{{\mathscr{L}}({{\boldsymbol{\xi }}}_{j+1})}{{\mathscr{L}}({{\boldsymbol{\xi }}}_{j})}=\exp \,\left(-\displaystyle \frac{{\chi }_{{{\boldsymbol{\xi }}}_{j+1}}^{2}-{\chi }_{{{\boldsymbol{\xi }}}_{j}}^{2}}{2}\right).\end{eqnarray} \tag{ 20 }$$

ExoSOFT also provides ${\mathscr{L}}$ and p in a format that easily interfaces with MCMC packages like emcee (Foreman-Mackey et al. 2013), using the combined function ${ln}({\mathscr{L}}({\boldsymbol{\xi }}))+{ln}(p({\boldsymbol{\xi }}))$. As such, emcee has been made a dependency of ExoSOFT and can be selected for use with a flag in the settings file.

We chose prior probabilities, or "priors," for each parameter with simplicity in mind. These should only be considered as suggestions, and may be easily changed in the ExoSOFT settings files to suit the user's preferences or to take into account the results of more recent work. The complete set of parameters in our full joint analysis are (m₁, m₂, ϖ, P, i, e, T_o, Ω, ω, γ). Independent of their distribution, the priors are all given wide boundary limits, each available for change within the settings. Those of Ω and ω are drawn without limits, but the stored values are then shifted into [0, 2π].

For m₁, the prior can be represented by the initial mass function (IMF) or present day mass function (PDMF). Chabrier (2003) found functional forms of these for a range of stellar populations, including those of the galactic disk. For m₂, however, particularly for stellar binaries, the IMF may be a poor approximation to the companion mass function (CMF); that of Metchev & Hillenbrand (2009) might be more appropriate. More recent investigations (e.g., Brandt et al. 2014) were consistent with Metchev & Hillenbrand (2009) down to a few Jupiter masses. PDMF and IMF can be found in Table 1 of Chabrier (2003), and that of the CMF in Equation (8) of Metchev & Hillenbrand (2009). By default we use the PDMF for the primary and CMF for the companion, with the IMF being an option also available to the user for either body.

The prior for the parallax (ϖ) may be derived from the assumption that stars are uniformly distributed in space. For nearby objects, within a few 100 pc, this prior is

$$\begin{eqnarray}&&p(\varpi )\propto \displaystyle \frac{1}{{\varpi }^{4}}.\end{eqnarray} \tag{ 21 }$$

In cases where high-S/N distance measurements are available, this information can be included with a Gaussian prior with center and width equal to its distance and error, respectively. These two components would then be multiplied to form a combined prior for the parallax. For distant objects, over a few hundred parsecs, the assumption of uniform stellar distribution is invalid. Low-S/N measurements of the parallax, in combination with Equation (21), would lead to unrealistically large distances. ExoSOFT currently requires accurately measured parallaxes. In practice, these are available for nearby bright stars from Hipparcos (van Leeuwen 2007), and the sample of objects with well-measured parallaxes will soon expand dramatically, thanks to Gaia.

The adopted prior for the period is a power law, $p(P)\propto {P}^{\gamma }$. Radial velocity surveys favor γ ≈ −0.7 (Cumming et al. 2008), while γ = −1 gives equal numbers per logarithmic period.

We assume systems to be randomly oriented with $p(i)\,\propto \sin (i)$.

Following the results of investigations into the eccentricity distributions of low-mass companions (e.g., Duquennoy & Mayor 1991; Cumming et al. 2008; Jurić & Tremaine 2008; Shen & Turner 2008; Hogg et al. 2010; Wang & Ford 2011; Kipping 2013; Kipping & Sandford 2016; Shabram et al. 2016), ExoSOFT offers multiple priors for the eccentricity. Currently, these include uniform, Rayleigh, Rayleigh+Exponential, Beta, $p(e)=2e$ (for companions with P > 1000 days), and those of Shen & Turner (2008). Additional priors can be added and modified manually by the user. This can be accomplished by placing a customized version of the priors module in the user's current working directory.

The remainder of the parameters, (T₀, Ω, ω, γ), has been given uniform priors.

Given an infinitely long Markov chain, convergence to the posterior probability distributions is guaranteed. In practice, when trying to integrate Bayes' theorem for a particular model, one of the primary goals is to minimize the number of samples necessary to achieve sufficient convergence. Therefore, before performing MCMC, a metric to measure its convergence is required along with steps to optimize the convergence rate.

It is difficult to determine if a Markov Chain has converged to the posterior probability distribution. From a given set of parameter values, there is no guarantee that the likelihood function will smoothly increase toward the global maximum of the posterior. In cases with complicated likelihood topography, there might be many places where a chain could become stuck. Chains started near a local maximum could falsely appear to have converged or take a long time to diffuse into other regions. For this reason, the specific starting position of a Markov chain, and any following samples drawn far from the likelihood peaks, must be removed to avoid any disproportionalities in the posterior distributions they might induce. This initial period is referred to as the "burn-in." Although the exact number of steps involved is not clearly defined, Tegmark et al. (2004) took it to be over when the likelihood in a particular chain was equal to the median likelihood of all chains combined.

The rate at which a Markov chain explores the posteriors is related to the width, or σ, of each parameter's proposal distribution . For large σ, proposed steps tend to lie far from the current one and in a region where they are unlikely to be accepted. By contrast, a small σ would lead to a high rate of acceptance, but a slow rate at which the posteriors are explored. Gelman et al. (1996, p. 599) showed that the highest diffusion speed of a one-dimensional model was achieved when the probability of accepting a proposed step was 0.441. They tested a variety of multi-dimensional models, concluding that an acceptance rate of ∼25% is suitable for higher dimensional models and ∼50% for those with only one or two dimensions.

With a Markov sampling scheme, the dependence between neighboring samples can dramatically reduce the number of effectively independent steps taken. Unfortunately, prior to running a Markov chain, there are no direct methods to predict how many independent steps will be taken over a certain number of samples. Instead, following the completion of all Markov chains, a common practice is to estimate if enough steps were taken using either the Gelman–Rubin statistic (Gelman & Rubin 1992) or the integrated autocorrelation time (Goodman & Weare 2010). As discussed in Foreman-Mackey et al. (2013) and Goodman & Weare (2010), the integrated autocorrelation time can be represented by

$$\begin{eqnarray}&&\tau =1+2\displaystyle \sum _{T=1}^{\infty }\,\displaystyle \frac{C(T)}{C(0)},\end{eqnarray} \tag{ 22 }$$

where T is the time lag necessary for covariance between samples, C(T), to go to zero, ensuring their independence. With this, the number of effectively independent steps in a Markov chain of length L_c is ∼L_c/τ.

ExoSOFT offers four of its own algorithms for exploring the parameter space, each for a different task: standard "shotgun" Monte Carlo, Simulated Annealing, Sigma Tuning, and MCMC. It is also ready to be interfaced with third-party sampling tools.

As discussed in Section 3.2, standard Monte Carlo is suitable for models with few free parameters, or when the likelihood function is not strongly peaked. In such cases, it could in fact sample the posterior probability distributions significantly faster than MCMC. With each sample being completely independent, it can also be useful for getting an idea of where the regions of minimum χ² lie.

In cases where standard Monte Carlo is inefficient, Markov sampling techniques are suggested. To minimize the burn-in, suitable starting positions for the Markov chains are found using Simulated Annealing. As discussed in Kirkpatrick et al. (1983), Simulated Annealing also makes use of the Metropolis–Hastings algorithm with an additional "temperature" factor in the likelihoods ratio, modifying this ratio in the rejection function to

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{\mathscr{L}}({{\boldsymbol{\xi }}}_{j+1})}{{\mathscr{L}}({{\boldsymbol{\xi }}}_{j})}\right)}^{1/T}=\exp \,\left(-\displaystyle \frac{{\chi }_{{{\boldsymbol{\xi }}}_{j+1}}^{2}-{\chi }_{{{\boldsymbol{\xi }}}_{j}}^{2}}{2T}\right).\end{eqnarray} \tag{ 23 }$$

At high temperatures, Equation (23) is temperature dominated, and likely to accept even unfavorable steps, thereby mitigating any local topographical bumps encountered. A cooling scheme started at a sufficiently high temperature (⪆100) explores a large portion of the parameter space and is slowly "annealed" toward the parameter set with the global maximum likelihood. In our implementation, the chains are repeatedly heated and cooled, each time from a lower starting temperature, while maintaining a list of the samples with the highest likelihoods achieved so far. We stop this process when the fits have all converged to a single set of parameter peaks in the posterior probability distributions.

After finding an initial set of parameters, ExoSOFT can sample the posterior probability distributions with either MCMC or the affine-invariant ensemble sampler built into emcee (Foreman-Mackey et al. 2013). When using MCMC, we first apply a version of the Metropolis–Hastings algorithm to determine the appropriate widths for each proposal distribution. Periodically, the acceptance rates are calculated, typically after ∼200 samples per varied parameter, and the fixed σ values of the normal proposal distributions for each parameter incremented accordingly until stable acceptance rates between 25%–35% are reached; this process is occasionally referred to as Sigma Tuning. In contrast to the earlier Simulated Annealing and later MCMC stages, Sigma Tuning only requires a single chain and has been implemented in such a way that it automatically completes once the acceptance rates stabilize. In cases where the observational data is sparse, or parameters highly correlated, the use of an ensemble sampler has been proven to increase the sampling efficiency (Goodman & Weare 2010; Foreman-Mackey et al. 2013). Thus, our recommended approach is to utilize the starting position produced by Simulated Annealing to initialize the walkers of the affine-invariant ensemble sampler built into emcee (Foreman-Mackey et al. 2013).

The solutions for determining a set of initial parameters and tuning the proposal distribution widths lend themselves well to integration with other third-party sampling tools, including PyMC (Fonnesbeck et al. 2015) and PyMultiNest (Feroz et al. 2009).

Our model was verified with synthetic data produced using an independent Keplerian tool built upon the PyAstronomy package.¹¹ The observable values had their errors realized from Gaussian distributions, using a percentage of each observable's absolute mean as the distribution widths. Such synthetic data were produced for various systems spanning the range of exoplanets and binary star systems observed to date to ensure ExoSOFT's ability to fit any system of interest.

For demonstration purposes, the approximate orbital elements of Jupiter were used assuming a host star at a distance of 20 pc, or a parallax of 50 mas, with the orbital plane inclined by 45° to Earth's line of sight. A variety of measurement errors were assumed for the same orbital elements to show the convergence of the resulting posterior distributions with improving data quality; the errors indicated are a percentage of each observable parameter's absolute mean value.

Figure 1 shows the marginalized posteriors using realizations of 10%, 5%, and 1% Gaussian noise on the measured parallax, radial velocity, and angular separation. The 68% confidence regions shrink by a factor of ∼10 from the 10% to 1% cases, shown as the green and red lines in Figure 1, respectively. The posterior distributions of T₀ and ω (not shown) have larger widths and decrease faster with improving data quality; this is an artifact of our choice of a nearly circular orbit. In all of the hypothetical Jupiter demonstration cases, ExoSOFT used a uniform prior for e, appropriate for planets. This was implemented by drawing ω and e simultaneously according to the parametrization from Albrecht et al. (2012).

Zoom In Zoom Out Reset image size

Figure 2 shows the resulting marginalized posteriors for one realization of 5% measurement errors; Table 2 provides a complete summary of the fit values. The covariance between a subset of the orbital elements is represented in a corner plot in Figure 4. Here the relationship between the masses and semimajor axis can be seen, arising from Kepler's third law in Equation (11) being directly included in the model equations.

Zoom In Zoom Out Reset image size

The best sample, defined as the set of stored parameters with the highest likelihood, has χ² = 46.9 for 35 degrees of freedom (25 epochs of radial velocity data, 10 epochs of astrometry in R.A. and decl., and 10 fitted parameters), for a reduced χ² = 1.34. Figure 3 shows the orbit for the best sample plotted atop the noisy astrometry and radial velocity data. The residuals, as expected, scatter about the best sample's orbit and typically lie ∼1σ away.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The joint model given in Section 2 mutually constrains the parameters P, m_tot, ω*, i, e, and T₀, where ω* represents the value for either body following the relation ω₁ = ω₂ + π. The other parameters (γ, m₂, Ω₂, ϖ) are constrained by a single observation type, or in the form of a combined parameter, as is the case for m_tot = m₁ + m₂. We performed the final fitting with the affine-invariant ensemble sampler of emcee and investigated the mixing and convergence of each parameter using the integrated autocorrelation time, discussed in Section 3.3 and represented in Equation (22). The least convergent parameter in all three error realizations was T₀, with a total of 5 × 10⁶ samples distributed over 1000 walkers, its τ was 690 for the 5% case. This slow convergence is an artifact of our choice of a nearly circular orbit (the definition of T₀ is arbitrary for a circular orbit) and is much less pronounced for more eccentric systems.

Although this test case is at the edge of current instruments' sensitivity, the ∼02 separation is within the expected limits of upcoming instruments like CHARIS (Peters et al. 2012), SPHERE (Claudi et al. 2006), and GPI (Macintosh et al. 2014). However, even with the improved contrast ratios expected at these separations, this hypothetical planetary system would need to be dramatically younger than Jupiter for detection through direct imaging techniques. The uncertainties in m₁ and m₂ also depend on the precision of the measured parallax, though this is already better than 1% for many bright, nearby stars (van Leeuwen 2007) and will further improve with the forthcoming Gaia data release.

An example of ExoSOFT applied to real data can be found in Hełminiak et al. (2016). In that paper, we introduce the first direct imaging detections of the companion to V450 And, solve for its joint orbital solution, and compare the resulting dynamical masses to those from stellar model fits based on the photometry.