Abstract
Markov Chain Monte Carlo (MCMC) methods for sampling probability density functions (combined with abundant computational resources) have transformed the sciences, especially in performing ...probabilistic inferences, or fitting models to data. In this primarily pedagogical contribution, we give a brief overview of the most basic MCMC method and some practical advice for the use of MCMC in real inference problems. We give advice on method choice, tuning for performance, methods for initialization, tests of convergence, troubleshooting, and use of the chain output to produce or report parameter estimates with associated uncertainties. We argue that autocorrelation time is the most important test for convergence, as it directly connects to the uncertainty on the sampling estimate of any quantity of interest. We emphasize that sampling is a method for doing integrals; this guides our thinking about how MCMC output is best used.
We derive analytic, closed form, numerically stable solutions for the total flux received from a spherical planet, moon, or star during an occultation if the specific intensity map of the body is ...expressed as a sum of spherical harmonics. Our expressions are valid to arbitrary degree and may be computed recursively for speed. The formalism we develop here applies to the computation of stellar transit light curves, planetary secondary eclipse light curves, and planet-planet/planet-moon occultation light curves, as well as thermal (rotational) phase curves. In this paper, we also introduce starry, an open-source package written in C++ and wrapped in Python that computes these light curves. The algorithm in starry is six orders of magnitude faster than direct numerical integration and several orders of magnitude more precise. starry also computes analytic derivatives of the light curves with respect to all input parameters for use in gradient-based optimization and inference, such as Hamiltonian Monte Carlo (HMC), allowing users to quickly and efficiently fit observed light curves to infer properties of a celestial body's surface map. (Please see https://github.com/rodluger/starry, https://rodluger.github.io/starry/, and https://doi.org/10.5281/zenodo.1312286).
The past two decades have seen a major expansion in the availability, size, and precision of time-domain data sets in astronomy. Owing to their unique combination of flexibility, mathematical ...simplicity, and comparative robustness, Gaussian processes (GPs) have emerged recently as the solution of choice to model stochastic signals in such data sets. In this review, we provide a brief introduction to the emergence of GPs in astronomy, present the underlying mathematical theory, and give practical advice considering the key modeling choices involved in GP regression. We then review applications of GPs to time-domain data sets in the astrophysical literature so far, from exoplanets to active galactic nuclei, showcasing the power and flexibility of the method. We provide worked examples using simulated data, with links to the source code; discuss the problem of computational cost and scalability; and give a snapshot of the current ecosystem of open-source GP software packages. In summary:
GP regression is a conceptually simple but statistically principled and powerful tool for the analysis of astronomical time series.
It is already widely used in some subfields, such as exoplanets, and gaining traction in many others, such as optical transients.
Driven by further algorithmic and conceptual advances, we expect that GPs will continue to be an important tool for robust and interpretable time-domain astronomy for many years to come.
We present an update to the EVEREST K2 pipeline that addresses various limitations in the previous version and improves the photometric precision of the light curves. We develop a fast regularization ...scheme for pixel-level decorrelation (PLD) and adapt the algorithm to include the PLD vectors of other stars to enhance the predictive power of the model and minimize overfitting, particularly for faint stars. We also modify PLD to work for saturated stars and improve its performance on variable stars, although some high-frequency variables may still suffer from overfitting. On average, EVEREST 2.0 light curves have 10-20% higher photometric precision than those in version 1, yielding the highest-precision light curves at all magnitudes of any publicly available K2 catalog. For most K2 campaigns, we recover the original Kepler precision to at least = 14, and to at least = 15 for campaigns 1, 5, 6, and 13. We also detrend most short-cadence targets observed by K2, obtaining even higher photometric precision for these stars. Like all aggressive, flexible models, EVEREST is prone to overfitting, and may cause a decrease in transit depths by ∼10%; we urge users to mask signals of interest using our open-source software, which we show removes this bias. Light curves for campaigns 0-8 and 10-13 are available online in the EVEREST catalog, which will be updated with future campaigns. EVEREST 2.0 is open source and is coded in a framework that can be adapted to other photometric surveys, including Kepler and the upcoming TESS mission.
Abstract
Variability in the light curves of spotted, rotating stars is often non-sinusoidal and quasi-periodic – spots move on the stellar surface and have finite lifetimes, causing stellar flux ...variations to slowly shift in phase. A strictly periodic sinusoid therefore cannot accurately model a rotationally modulated stellar light curve. Physical models of stellar surfaces have many drawbacks preventing effective inference, such as highly degenerate or high-dimensional parameter spaces. In this work, we test an appropriate effective model: a Gaussian Process with a quasi-periodic covariance kernel function. This highly flexible model allows sampling of the posterior probability density function of the periodic parameter, marginalizing over the other kernel hyperparameters using a Markov Chain Monte Carlo approach. To test the effectiveness of this method, we infer rotation periods from 333 simulated stellar light curves, demonstrating that the Gaussian process method produces periods that are more accurate than both a sine-fitting periodogram and an autocorrelation function method. We also demonstrate that it works well on real data, by inferring rotation periods for 275 Kepler stars with previously measured periods. We provide a table of rotation periods for these and many more, altogether 1102 Kepler objects of interest, and their posterior probability density function samples. Because this method delivers posterior probability density functions, it will enable hierarchical studies involving stellar rotation, particularly those involving population modelling, such as inferring stellar ages, obliquities in exoplanet systems, or characterizing star–planet interactions. The code used to implement this method is available online.
emcee: The MCMC Hammer Foreman-Mackey, Daniel; Hogg, David W.; Lang, Dustin ...
Publications of the Astronomical Society of the Pacific,
03/2013, Letnik:
125, Številka:
925
Journal Article
Recenzirano
Odprti dostop
ABSTRACT We introduce a stable, well tested Python implementation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) proposed by Goodman & Weare (2010). The code is open ...source and has already been used in several published projects in the astrophysics literature. The algorithm behind emcee has several advantages over traditional MCMC sampling methods and it has excellent performance as measured by the autocorrelation time (or function calls per independent sample). One major advantage of the algorithm is that it requires hand-tuning of only 1 or 2 parameters compared to ∼N2 for a traditional algorithm in an N-dimensional parameter space. In this document, we describe the algorithm and the details of our implementation. Exploiting the parallelism of the ensemble method, emcee permits any user to take advantage of multiple CPU cores without extra effort. The code is available online at http://dan.iel.fm/emcee under the GNU General Public License v2.
No true extrasolar Earth analog is known. Hundreds of planets have been found around Sun-like stars that are either Earth-sized but on shorter periods, or else on year-long orbits but somewhat ...larger. Under strong assumptions, exoplanet catalogs have been used to make an extrapolated estimate of the rate at which Sun-like stars host Earth analogs. These studies are complicated by the fact that every catalog is censored by non-trivial selection effects and detection efficiencies, and every property (period, radius, etc.) is measured noisily. Here we present a general hierarchical probabilistic framework for making justified inferences about the population of exoplanets, taking into account survey completeness and, for the first time, observational uncertainties. We are able to make fewer assumptions about the distribution than previous studies; we only require that the occurrence rate density be a smooth function of period and radius (employing a Gaussian process). By applying our method to synthetic catalogs, we demonstrate that it produces more accurate estimates of the whole population than standard procedures based on weighting by inverse detection efficiency. We apply the method to an existing catalog of small planet candidates around G dwarf stars. We confirm a previous result that the radius distribution changes slope near Earth's radius. We find that the rate density of Earth analogs is about 0.02 (per star per natural logarithmic bin in period and radius) with large uncertainty. This number is much smaller than previous estimates made with the same data but stronger assumptions.
Abstract
Stellar obliquity, the angle between a planet’s orbital axis and its host star’s spin axis, traces the formation and evolution of a planetary system. In transiting-exoplanet observations, ...only the sky-projected stellar obliquity can be measured, but this can be deprojected using an estimate of the stellar obliquity. In this paper, we introduce a flexible, hierarchical Bayesian framework that can be used to infer the stellar obliquity distribution solely from sky-projected stellar obliquities, including stellar inclination measurements when available. We demonstrate that while a constraint on the stellar inclination is crucial for measuring the obliquity of an individual system, it is not required for robust determination of the population-level stellar obliquity distribution. In practice, the constraints on the stellar obliquity distribution are mainly driven by the sky-projected stellar obliquities. When applying the framework to all systems with measured sky-projected stellar obliquity, which are mostly hot Jupiter systems, we find that the inferred population-level obliquity distribution is unimodal and peaked at zero degrees. Misaligned systems have nearly isotropic stellar obliquities with no strong clustering near 90°. The diverse range of stellar obliquities prefers dynamic mechanisms, such as planet–planet scattering after a convergent disk migration, which could produce both prograde and retrograde orbits of close-in planets with no strong inclination concentrations other than that at 0°.
Fast Direct Methods for Gaussian Processes Ambikasaran, Sivaram; Foreman-Mackey, Daniel; Greengard, Leslie ...
IEEE transactions on pattern analysis and machine intelligence,
2016-Feb.-1, 2016-Feb, 2016-2-1, 20160201, Letnik:
38, Številka:
2
Journal Article
Recenzirano
Odprti dostop
A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and ...variance simply requires the evaluation of the corresponding Gaussian density. In the n-dimensional setting, however, it requires the inversion of an n x n covariance matrix, C, as well as the evaluation of its determinant, det(C). In many cases, such as regression using Gaussian processes, the covariance matrix is of the form C = σ 2 I + K, where K is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix C is typically dense, causing standard direct methods for inversion and determinant evaluation to require O(n 3 ) work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix C can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an O(n log 2 n) algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant det(C), permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining K. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.
Among the available methods for dating stars, gyrochronology is a powerful one because it requires knowledge of only the star's mass and rotation period. Gyrochronology relations have previously been ...calibrated using young clusters, with the Sun providing the only age dependence, and are therefore poorly calibrated at late ages. We used rotation period measurements of 310 Kepler stars with asteroseismic ages, 50 stars from the Hyades and Coma Berenices clusters and 6 field stars (including the Sun) with precise age measurements to calibrate the gyrochronology relation, whilst fully accounting for measurement uncertainties in all observable quantities. We calibrated a relation of the form P = A
n
× (B − V − c)
b
, where P is rotation period in days, A is age in Myr, B and V are magnitudes and a, b and n are the free parameters of our model. We found
$a = 0.40^{+0.3}_{-0.05}$
,
$b = 0.31^{+0.05}_{-0.02}$
and
$n = 0.55^{+0.02}_{-0.09}$
. Markov Chain Monte Carlo methods were used to explore the posterior probability distribution functions of the gyrochronology parameters and we carefully checked the effects of leaving out parts of our sample, leading us to find that no single relation between rotation period, colour and age can adequately describe all the subsets of our data. The Kepler asteroseismic stars, cluster stars and local field stars cannot all be described by the same gyrochronology relation. The Kepler asteroseismic stars may be subject to observational biases; however, the clusters show unexpected deviations from the predicted behaviour, providing concerns for the overall reliability of gyrochronology as a dating method.
PDF
