Improving \textsl{Gaia} parallax precision with a data-driven model of stars

TL;DR

This paper introduces a data-driven, empirical prior for stellar distances using Gaia and 2MASS data, significantly improving parallax precision by leveraging a large, noise-deconvolved color-magnitude model.

Contribution

It develops a novel noise-deconvolved empirical prior based on 2MASS colors and Gaia data, enhancing parallax estimates without relying on stellar physics or Galactic models.

Findings

Median parallax precision improved by 20% over Gaia TGAS

14% of estimates are more than twice as precise

Validated with star cluster M67, revealing clearer cluster structure

Abstract

Converting a noisy parallax measurement into a posterior belief over distance requires inference with a prior. Usually this prior represents beliefs about the stellar density distribution of the Milky Way. However, multi-band photometry exists for a large fraction of the \textsl{\small{Gaia}} \textsl{\small{TGAS}} Catalog and is incredibly informative about stellar distances. Here we use \textsl{\small{2MASS}} colors for 1.4 million \textsl{\small{TGAS}} stars to build a noise-deconvolved empirical prior distribution for stars in color--magnitude space. This model contains no knowledge of stellar astrophysics or the Milky Way, but is precise because it accurately generates a large number of noisy parallax measurements under an assumption of stationarity; that is, it is capable of combining the information from many stars. We use the Extreme Deconvolution (\textsl{\small{XD}}) algorithm---an Empirical Bayes approximation to a full hierarchical model of the true parallax and photometry of every star---to construct this prior. The prior is combined with a \textsl{\small{TGAS}} likelihood to infer a precise photometric parallax estimate and uncertainty (and full posterior) for every star. Our parallax estimates are more precise than the \textsl{\small{TGAS}} catalog entries by a median factor of 1.2 (14% are more precise by a factor >2) and are more precise than previous Bayesian distance estimates that use spatial priors. We validate our parallax inferences using members of the Milky Way star cluster M67, which is not visible as a cluster in the \textsl{\small{TGAS}} parallax estimates, but appears as a cluster in our posterior parallax estimates. Our results, including a parallax posterior pdf for each of 1.4 million \textsl{\small{TGAS}} stars, are available in companion electronic tables.

Figures

Equations

$$y_{\mathrm{true},n}$$

$$y_n$$

$$\displaystyle p(y_{n}\,|\,m,b,t^{2},\sigma_{n}^{2})\$$

$$=$$

$$p(\{y_n\}|m,b,t^2,{\sigma }_n^2)=\prod _{n}p(y_n|m,b,t^2,{\sigma }_n^2).$$

$${\chi }^2=\sum _n\frac{[y_n-[m{x}_n+b]{]}^2}{[{\sigma }_n^2+t_n^2]}+\sum _n\ln (2\pi [{\sigma }_n^2+t^2]).$$

$$p(y_{\mathrm{true},n}|y_n,{\tilde{\sigma }}_n)$$

$$\mathrm{where}p(y_n|y_{\mathrm{true},n},{\sigma }_n)$$

$$\mathrm{and}p(y_{\mathrm{true},n})$$

$$p({\varpi }_{\mathrm{true}}|\varpi ,{\sigma }_{\varpi }^2)=\frac{1}{Z}p(\varpi |{\varpi }_{\mathrm{true}},{\sigma }_{\varpi }^2)p({\varpi }_{\mathrm{true}}),$$

$$y_0=\varpi 10^{\frac{1}{5}{J}_n^{\mathrm{c}}}=10^{\frac{1}{5}{M}_{J,n}^{\mathrm{c}}+2}$$

$$\mathrm{true}\mathrm{vector}$$

$$\mathrm{data}\mathrm{vector}$$

$$\mathrm{with}\mathrm{apparent}\mathrm{magnitude}$$

$$\mathrm{and}\mathrm{color}$$

$$\mathrm{under}\mathrm{the}\mathrm{assumptions}$$

$$p({\mathbf{y}}_n|{\mathbf{Y}}_n,I_n)=\mathcal{N}({\mathbf{y}}_n|{\mathbf{Y}}_n,{\mathbf{C}}_n)$$

$${\mathbf{C}}_n=\left[\begin{array}{cc}{\left({\sigma }_{\varpi }10^{0.2J_n^{\mathrm{c}}}\right)}^2& 0\\ 0& {\sigma }_{J,n}^2+{\sigma }_{K,n}^2\end{array}\right]$$

$$p({\mathbf{Y}}_n|{\mathbf{y}}_n,I_n)\propto p({\mathbf{y}}_n|{\mathbf{Y}}_n,I_n)p({\mathbf{Y}}_n)$$

$$p\bigl{(}\ \boldsymbol{{Y}}\,|\,\left\{A_{k},\boldsymbol{{\mu}}_{k},{\mathbf{{V}}}_{k}\right\}_{k=1}^{K}\bigr{)}=\sum_{k=1}^{K}A_{k}\ \mathcal{N}(\boldsymbol{{Y}}\,|\,\boldsymbol{{\mu}}_{k},{\mathbf{{V}}}_{k})$$

$$\hat{\alpha }$$

$$p({\varpi }_{\mathrm{true}}|\mathbf{y},\hat{\alpha })$$

$$\mathrm{where}f({\varpi }_{\mathrm{true}})$$

Full text

catalogs — Hertzsprung–Russell and C–M diagrams — methods: statistical — parallaxes — stars: distances

††software: The code used in this project is available from https://github.com/andersdot/photoParallax. This research utilized the following open-source Python packages: Astropy (Astropy Collaboration et al. 2013), matplotlib (Hunter 2007), and numpy (Van der Walt et al. 2011). This work additionally used the Gaia science archive (https://gea.esac.esa.int/archive/).

1 Introduction

The Gaia Mission (Gaia Collaboration et al. 2016b) will soon deliver more than a billion stellar parallax measurements. Only a small fraction (but large number) of these measurements will be precise and purely astrometric. Gaia will use astrometric parallaxes to determine the distances of the more precise stars, and then calibrate spectrophotometric models (Bailer-Jones et al. 2013). These spectrophotometric models, along with Gaia’s on-board low-resolution $BP/RP$ spectrophotometry, can be used to provide more precise parallax estimates for stars with poor astrometric parallax data. The full stack required for these parallax inferences is complex. It involves modeling the stars, as well as the dust, in the Milky Way, and the response of the telescope itself. It might be possible to eliminate some of this complexity using data driven models which don’t rely on physical models.

Projects like The Cannon (Ness et al. 2015; Casey et al. 2016; Ho et al. 2017) and Avast (Bedell et al. in preparation), explore the extent to which our predictive models of stars could be purely data-driven or statistical. That is, under what circumstances could the data themselves deliver more precise or more accurate information than any theoretical or physics-based model? The answer to this question is extremely context-dependent: It depends on what data are available, how much data are available, and what questions are being asked. In general, data-driven models contain fewer assumptions than physics-driven models, but they are also usually less interpretable; the model presented here has these same properties. Yet, the Gaia data set is ideal for thinking about these kinds of purely statistical models of stars. The different stars observed by Gaia are measured at very different signal-to-noise ratios, and a data-driven model of the stellar color–magnitude distribution can capitalize on the precise information coming from the high signal-to-noise stars and deliver valuable information about the low signal-to-noise stars. This information can then deliver precise parallax (or distance or absolute magnitude) inferences for all of the stars.

Here we present a demonstration-of-concept that shows that a data-driven model of the Gaia data could in principle deliver very precise parallax measurements, without any input of the physics or numerical models of stellar evolution, interiors, or photospheres. That is, this is a data-driven or purely empirical photometric parallax method. A photometric parallax is a parallax (or distance) estimate based on the photometry (colors and apparent magnitude or magnitudes) of a star. This idea was pioneered by Jurić et al. (2008) using the Sloan Digital Sky Survey (York et al. 2000). Fundamentally, photometric parallaxes capitalize on the fact that the apparent magnitude of a star is a strong function of distance, and that the absolute magnitude is a strong function of color (or surface temperature) and evolutionary stage.

Most stellar models used to generate photometric parallaxes have been physical models, based on gravity, fluid mechanics, radiative transfer, nuclear reactions, and atomic and molecular transitions. Because there are small issues with each of these components, the physical models are inaccurate in detail, and therefore produce parallax (and distance) estimates that are biased. In addition, they build a long list of physical assumptions (about nuclear reactions, convection, and thermal timescales, for example) into the parallax estimates. When using physical models, it is impossible to deliver photometric parallax estimates that involve minimal assumptions.

However, the use of physical models for distance estimation is not necessary. It is possible to build a stellar photometric model from the data themselves, because there are stars at a wide range of luminosities with useful parallax information, and this range of stellar luminosities with useful parallax information will grow with future Gaia data releases. One challenge is that different stars are observed at different levels of parallax precision, so it requires relatively sophisticated technology to build this data-driven model using all of the data available, fairly and responsibly.

Here we build and use just such a data-driven model to infer more precise parallaxes for each star. In particular, we use all of the Gaia TGAS data (Gaia Collaboration et al. 2016a; these data were generated by the methodology of Michalik et al. 2015), that match to the 2MASS photometry (Skrutskie et al. 2006) and lie in the footprint of PanSTARRS1 (PS1; Chambers et al. 2016) to make use of the Green et al. (2015) dust map. We use this data to make a model of the noise-free (or low-noise) color–magnitude diagram (CMD) of stars. We build the CMD using an Empirical-Bayes approach known as Extreme Deconvolution (XD; Bovy et al. 2011a). This method deconvolves heteroskedastic data to derive an estimate of the noise-free or low-noise distribution that would have been observed with far better data. This method has been used in astrophysics previously to model the velocity distribution of stars in the Solar Neighborhood (Hogg et al. 2005; Bovy et al. 2009) and to perform quasar target selection in Sloan Digital Sky Survey data (Bovy et al. 2011b; Bovy et al. 2012). Its advantage over other deconvolution methods is that it takes as input, and handles in a principled way, heteroskedastic data. Its disadvantage, for the present purposes, is that it requires or implies a strictly Gaussian noise model. As we show below, this requires us to transform the CMD to a space in which the noise is approximately Gaussian.

We then use the XD output—the deconvolved CMD—as a prior for use in inference of individual stellar parallaxes. These inferences, one per star, provide much more precise parallax, distance, or absolute magnitude estimates than we get from each star’s primary TGAS parallax measurement alone. That is, we are using the XD model to de-noise the Gaia parallax measurements. Technically, since XD produces a maximum-marginalized-likelihood estimator of the deconvolved distribution, its use as a prior is technically using the data twice and is therefore only approximately valid. However, since the data set is large, the XD approximation is not bad; it is sometimes known as the “empirical Bayes” method, and is well studied.111For a primer on “empirical Bayes”, see https://en.wikipedia.org/wiki/Empirical_Bayes_method.

A probabilistically valid approach that wouldn’t re-use the data is to perform a full hierarchical Bayesian inference. Exploration of that possibility, and its computational tractability, is among our long-term goals and motivations. Along those lines, this paper can be seen as a companion to a related project (Leistedt & Hogg, 2017), in which we use a slightly less appropriate model for the CMD, but in which we take a much more principled approach to the inference.

Although not a fully Bayesian hierarchical inference, the methodology presented in this paper makes exceedingly weak assumptions about the Milky Way, and similarly weak assumptions about the properties of stars. To date, there has been some discussion in the literature about how to use a parallax responsibly to infer a distance (Astraatmadja & Bailer-Jones 2016a; Astraatmadja & Bailer-Jones 2016b). These papers use a Gaia measurement to construct a parallax likelihood, and combine this with a distance prior. While this is sensible and correct, all present methods proposed along these lines build in informative (and known-to-be-wrong) assumptions about the line-of-sight distance distribution to Gaia stars. That is, they build in strong and informative assumptions about the Milky Way. The strongest new assumption made in this work is that for every star in TGAS, there are other, photometrically (and bolometrically) similar stars. That is a much weaker assumption than has been made in most other Bayesian distance estimations.

This project is fundamentally a demonstration of concept: We are only using a subset of the “small” (relative to the expected final Gaia data set) TGAS Catalog. We are performing only an approximation to full Bayesian hierarchical inference. We have to use photometry that is ground-based, rather than the full-precision, space-based photometry that the Gaia Mission will deliver. However, as we show below, we get very good results in terms of parallax precision. It is promising that we will be able to obtain even better results using later Gaia data releases, and eventually infer distances for $>10^{9}$ Gaia stars using photometric parallax methodologies but without any commitment to—or even use of—physical models of stars or the Milky Way.

2 Why does this work?

In this paper we use the noise-deconvolved CMD as a prior to infer more precise parallax probability distribution functions (PDFs) from Gaia measurements. Understanding the validity of this method can be challenging because we are using the data both for the prior and the likelihood in our inference. To ease understanding, here we demonstrate the general methodology using a simpler toy model with simulated data. Instead of inferring parallax PDFs using the 2D distribution of the CMD as our prior, for this simpler model we will infer posterior PDFs of the value $y$. We perform this simpler inference using a prior that is also a 2D distribution: a 1D Gaussian in the $y$ direction, with a running mean $\mu=m\,x+b$ and some thickness $t$. Here the thickness of the toy distribution mimics the true thickness, or intrinsic scatter, of the CMD due to age, mass, and metallicity. The $1\sigma$ contours as well as (noise free) samples from the toy model distribution are shown in Figure 1, upper-left panel in red. Similarly, noisy samples from this distribution are shown in the upper-right panel in grey and black. Using a representation of the true distribution (which we learn from the noisy data) as our prior, we can infer some posterior belief over $y_{\mathrm{true},n}$, a true sample from the distribution, from $y_{n}$, a sample from the noisy distribution.

Setting up the problem in more detail, we have the true samples from the underlying distribution and the noisy samples

[TABLE]

where $\Delta_{n}$ is drawn from the Gaussian function $\mathcal{N}(\Delta_{n}\,|\,0,t^{2})$—which accounts for the intrinsic width—and $\epsilon_{n}$ is drawn from the Gaussian function $\mathcal{N}(\epsilon_{n}\,|\,0,\sigma_{n}^{2})$—the noise distribution—so $\sigma_{n}$ is the measurement uncertainty. These represent true and noisy draws from a distribution that is a straight line with some thickness $t$, and are visualized in the upper panels of Figure 1. The upper-left shows 1024 realizations of the true values, and the upper-right, as well as lower-left, show 1024 realizations of the noisy values $y_{n}$ and their associated measurement uncertainties $\sigma_{n}$ as the grey points; we highlight 5 random points in black.

To infer a posterior PDF of the true $y$ value for each measured $y$ value we need a prior. We use an estimate of the underlying true distribution that the noisy data $y_{n}$ are drawn from, which we build from the noisy $y_{n}$ data themselves. We learn the noise-deconvolved function using marginalized, maximum-likelihood estimation (MMLE), which marginalizes over the true data values and maximizes the likelihood of the data. Here our likelihood for each individual datum is

[TABLE]

where we have marginalized over the true sample $y_{\mathrm{true},n}$. The resulting likelihood is a Gaussian because the noise model and our assumptions about the intrinsic width are both Gaussian. The likelihood of the full data set $\{y_{n}\}$ is then the product of the individual likelihoods

[TABLE]

Maximizing this joint likelihood of all the data is the same as minimizing

[TABLE]

The MMLE gives us the best fit parameters $\hat{m}$, $\hat{b}$, and $\hat{t}$. In Figure 1, lower-left panel, the true functions $y_{\mathrm{true}}=m_{\mathrm{true}}\,x+b_{\mathrm{true}}\pm t_{\mathrm{true}}$ are shown as the red lines, and the best fit MMLE functions $\hat{y}_{\mathrm{true}}=\hat{m}\,x+\hat{b}\pm\hat{t}$ are shown as the green lines. Now that we have an estimate of the true underlying distribution $\hat{y}_{\mathrm{true}}$, we use this as a prior to infer a posterior belief of the true values for each datum $p(y_{\mathrm{true},n}\,|\,y_{n},\tilde{\sigma}_{n}^{2})$.

Here is a good point to step back and acknowledge that our method is a bit strange. We are using the data to determine the noise-deconvolved distribution (an estimate of the true distribution), and then using that noise-deconvolved distribution to make inferences about the data. We are technically using the data twice! To be fully, technically correct, we should either learn the true distribution while leaving one datum out, the one datum we want to infer the true y value for, or else perform a full hierarchical Bayesian inference, in which we treat the CMD parameters probabilistically along with everything else. With large data sets, the MMLE or Empirical Bayes approximation is safe, in some sense, because no individual data point makes a large difference to the inference.

Returning to the inference: We use Bayes’ theorem to turn the likelihood of the noisy data into a posterior belief over the true $y$ value:

[TABLE]

In Figure 1 lower-left panel, the blue points represent the posterior beliefs (visualized as the expectation value and $1\sigma$ uncertainties) of the true values $y_{\mathrm{true},n}$ for the 5 highlighted (black) points from the upper-right panel. The red points are the associated true $y$ values. In the lower-right panel of Figure 1, the blue points represent the posterior beliefs, again visualized as the expectation value and $1\sigma$ uncertainty, for all 1024 $y_{n}$. When the true values are inferred in this way, the associated variance $\tilde{\sigma}_{n}^{2}$ is the convolution of the true distribution variance, $t^{2}$, and the measurement variance, $\sigma_{n}^{2}$. The points are colored by their measurement variance $\sigma_{n}^{2}$, which shows that noisier measurements are more influenced by the prior and get pulled closer to the center of the MMLE distribution, our prior. In contrast, more certain measurements remain closer to their measured values, $y_{n}$. We now have better estimates of each true $y$ value by inferring each true $y$ value from the measured $y$ value using an estimate of the true distribution the measured $y$ value was drawn from as the prior.

In detail, the posterior points shown in the lower-right panel of Figure 1 look more concentrated (vertically) than the true values shown in the upper-left panel. This is because the points in the lower-right (posterior) panel show not typical values or typical draws from the posterior; they show the modes (and, with error bars, the variances) of the posterior PDFs. When a measured value has a large variance, the posterior ends up close to the prior, which has a mode at the center line; the posterior modes therefore concentrate more towards the center line than any representative sampling. A plot of posterior samplings (as opposed to a plot of posterior modes) would show a scatter like that in the upper-left (truth) panel.

Bringing the toy model back to our original problem: We can get better estimates of the true parallax for each star by inferring the true parallax from the measured parallax and using the noise-deconvolved CMD, an estimate of the true CMD, as our prior. For the remainder of the paper, we will also maintain the color scheme presented here (where possible): truth will be red, data will be black, prior information will be green, and posterior information will be blue.

3 Assumptions and methods

To infer the true parallax for a star using the noise-deconvolved CMD as a prior, we make the following assumptions. The method we present is correct and justifiable under these assumptions, each of which is questionable in its own right. We return to criticize these assumptions—and the method that flows from them—in the discussion section below.

stationarity

We make an assumption of stationarity; that is, that the stars measured at lower signal-to-noise have analogs measured at higher signal-to-noise. In detail, this assumption is wrong. For example, more distant stars will have lower signal to noise ratio measurements and will also be different than local disk stars in age and metallicity. However, the stationarity assumption is fairly weak: the requirement is that, locally in the CMD, there is support among the high signal-to-noise stars for the types of stars also seen at low signal-to-noise. If this requirement is met, the results will not be strongly biased. Although this is a weak assumption, it is deep and fundamental to this work.

selection

We have no model for the selection function of TGAS (but see Bovy (2017)), nor selection volumes for any kinds of stars. For this reason, we are building a model of the contents of TGAS, not the properties of the stars in the Milky Way (nor any volume-limited component or region thereof). For this reason, our CMD prior might not be appropriate to use with other surveys.

big data

We assume that we have large numbers of stars, large enough that empirical-Bayes, or maximum-marginalized-likelihood, is a safe approximation to full Bayesian inference. This assumption will be least true in the least-populated parts of the CMD. In the extreme case, if a portion of the CMD has only one (statistically distinct) exemplar—for example, a white dwarf—inference for these cases will be biased.

noise model

We assume that the TGAS parallax uncertainties dominate the noise in any estimates of absolute magnitude. We further assume that the parallax and color uncertainties are Gaussian in form with correctly known variances. In detail, we assume that the parallax uncertainty estimates, (the parallax_error entry in the TGAS catalog), are correct. Here we are making use of the fact that when de-noising data, under estimating measurement noise is conservative, and leads to under-deconvolution. In other words, the features of the inferred CMD will typically be broader and more conservatively estimated than if we over-estimated the noise.

dust

We treat the PS1-based three-dimensional dust maps (Green et al. 2015) as correct in their median dust estimates, and their effects on color and magnitude. As we optimize the empirical Bayes model, we iteratively update each star’s dust correction at its updated, inferred distance. Our assumption is that these corrections are good enough, and uncertainties in these are not dominant, for neither color nor absolute magnitude.

mixture of Gaussians

We assume that the CMD can be represented by a mixture of Gaussians in a particular transformed space. This mixture we fix (with only heuristic analysis) at 128 Gaussians.

no physics

We make no use of stellar models; our only assumptions about stellar physics are generic and implicit (for example, that the CMD is somehow smooth and compact).

3.1 Setting up the Inference

Under these assumptions, we would like to infer the true parallax, $\varpi_{\mathrm{true}}$, for each star given its measured parallax, $\varpi$, and associated uncertainty, $\sigma_{\varpi}$, from TGAS. A first application of Bayes theorem leads to

[TABLE]

where $p(\varpi\,|\,\varpi_{\mathrm{true}},\sigma^{2}_{\varpi})$ is the likelihood of the observed parallax, $p(\varpi_{\mathrm{true}})$ is a prior PDF for the true parallax, and $Z$ is a normalization constant (the model evidence, or marginalized likelihood). The strength of our inference lies in how we generate our prior. Instead of assuming a functional form for the spatial distribution of stars in the Milky Way (Astraatmadja & Bailer-Jones 2016b), or relying on stellar models (Gaia Collaboration et al. 2016b), we build a noise-deconvolved estimate of the CMD from the observed colors and absolute magnitudes of our full set of stars. This is because supplementary to the trigonometric measure of distance, the photometry of stars is also informative of their distances. The intrinsic temperatures and luminosities of stars tightly correlate, leading to the observed colors and absolute magnitudes of stars also tightly correlating, usually visualized as the CMD. The tight correlation of the CMD benefits our inference because the luminosity or absolute magnitude of a star $M$, when compared with its brightness or apparent magnitude $m$, contains distance information $d$ via the fundamental relation $m=M+5\log(d/10\,\mathrm{pc})$. We include this in the inference by using the CMD as part of our prior for $\varpi_{\mathrm{true}}$. Specifically, we use the CMD generated from the data themselves, as detailed in the next section.

Equation 10 does not fully specify this procedure. Because our prior comes from the CMD, our inference will include not only the parallax from TGAS, but also the photometry of the stars. Specifically, we use the $J$ and $K_{s}$ band photometry from 2MASS, as well as the parallax from TGAS. Our true, latent values are the $(J-K_{s})^{\mathrm{true}}_{n}$ color and the $J$ band absolute magnitude $M^{\mathrm{true}}_{J,n}$ for each star, which we encapsulate in the true vector $\boldsymbol{{Y}}$. The observed data is the dust corrected color $(J-K_{s})^{\mathrm{c}}_{n}$, the TGAS parallax $\varpi_{n}$ in units of $[mas]$, and the dust corrected $J$ band magnitude $J^{\mathrm{c}}_{n}$, which we encapsulate in the data vector $\boldsymbol{{y}}$.

We now describe a singularity of our method, which allows us to conservatively deconvolve the data into a CMD prior without requiring a full hierarchical inference framework (which was explored in Leistedt & Hogg 2017). Instead of transforming the parallax and apparent magnitude into an absolute magnitude, the data vector $\boldsymbol{{y}}$ contains a transformation of the dust corrected absolute magnitude

[TABLE]

where, for this expression, it is important to remember that $\varpi$ is assumed to be in units of $[mas]$. The transformation of the (dust corrected) absolute magnitude is required to keep the uncertainties Gaussian, and Gaussian uncertainties are required to use XD, see below. Our prior assumptions $I$ are the uncertainties in each of these measurements, as well as the 3D dust map. So our full terminology is

[TABLE]

where $E(B-V)_{n}$ is the output of the dust model of Green et al. (2015), and $Q_{\lambda}$ is the correction factor for a photometric band assuming some reddening law.

We assume that the parallax likelihood is Gaussian, as well as the $(J-K_{s})^{\mathrm{c}}_{n}$ color likelihood. So our full likelihood for the $n$th star is

[TABLE]

with

[TABLE]

where $\mathcal{N}(\boldsymbol{{\mu}},{\mathbf{{\Sigma}}})$ represents a normal distribution with mean vector $\boldsymbol{{\mu}}$ and covariance matrix ${\mathbf{{\Sigma}}}$. Here we have assumed that the parallax uncertainty is significantly larger than the $J$-band magnitude uncertainty and is therefore the dominant uncertainty in the transformed absolute magnitude. In detail, we take $J_{n}^{\mathrm{c}}$ as a point estimate and only propagate the parallax uncertainty.

We can now write the multidimensional posterior expression analogous to Equation 10,

[TABLE]

where $p(\boldsymbol{{y}}_{n}\,|\,\boldsymbol{{Y}}_{n},I_{n})$ is the likelihood of the data vector, and $p(\boldsymbol{{Y}}_{n})$ is our prior. Before we describe how this can be turned into a posterior belief on the true parallax, we turn our attention to the prior, which will be generated from the data themselves. For this, we use XD, which takes advantage of the Gaussian likelihood.

3.2 Generating the Prior

To generate the empirical prior $p(\boldsymbol{{Y}}_{n})$, we fit our data in the transformed CMD space (see previous sub-section) using XD, which deconvolves the data $\{\boldsymbol{{Y}}_{n}\}$ with heterogeneous, Gaussian noise variances $\{{\mathbf{{C}}}_{n}\}$. XD generates an estimate of the underlying distribution from which the uncertain data were drawn, modeled as a mixture of $K$ Gaussians:

[TABLE]

where $\boldsymbol{{\mu}}_{k}$ and ${\mathbf{{V}}}_{k}$ are a two dimensional vector and matrix, respectively, which correspond to the location and covariance of the $k$th component (a Gaussian in our transformed CMD space). The relative amplitudes of the components are given by the scalars $\{A_{k}\}$.

Given the data $\{\boldsymbol{{Y}}_{n},{\mathbf{{C}}}_{n}\}$, XD maximizes the marginalized likelihood of the data, marginalizing over the true data. Assuming $\alpha=\left\{A_{k},\boldsymbol{{\mu}}_{k},{\mathbf{{V}}}_{k}\right\}_{k=1}^{K}$, it finds the coefficients

[TABLE]

which is made analytically tractable with the Gaussian likelihood and mixture-of-Gaussians prior (Bovy et al., 2011a).

DWH: write something about how our prior is a prior on distance

3.3 Back to the Inference

To infer $\varpi_{\mathrm{true}}$ once the CMD prior is constructed, we first do the inference in the 2D space of the CMD prior $p(\boldsymbol{{Y}}\,|\,\alpha)$, using Equation 14. We then project the 2D posterior into the 1D $\varpi_{\mathrm{true}}$ space by multiply by the Jacobian to calculate the posterior PDF over $\varpi_{\mathrm{true}}$

[TABLE]

with $f(\varpi_{\mathrm{true}})$ being a window function to insure $\varpi_{\mathrm{true}}$ is positive and to put $\varpi_{\mathrm{true}}$ on a similar grid for all the posteriors.

4 Data and results

We use stars crossmatched in TGAS and 2MASS. The match was done using a nearest-neighbor algorithm with a search radius of $4$ arcsec222http://portal.nersc.gov/project/cosmo/temp/dstn/gaia/tgas-matched-2mass.fits.gz. We also required that the stars lie within the observing footprint of PS1 to access the Green et al. (2015) 3D dust model. We require that the photometry have real values, and nonzero, real errors, and remove a small selection of 2MASS stars that have zero color and zero $J$-band magnitude. The full data set is visualized in the CMD in Figure 2. The left panel shows the point estimates of the color and absolute magnitude, using the point estimate of the parallax from TGAS. The $1\sigma$ and $2\sigma$ contours of the distribution are shown in grey. The right panel shows a subset of the data with the associated error bar for each star. The uncertainties in the colors are fairly well behaved, but the large uncertainties in some absolute magnitudes are due to the large uncertainties in the parallax measurements.

4.1 Dust

To generate the prior and evaluate the likelihood, we need to correct the 2MASS photometry for dust extinction. With the emergence of 3D dust maps, it is now possible to estimate dust corrections to stars within the Milky Way. The challenge is getting a proper estimate of the distance to the star before inferring $\varpi_{\mathrm{true}}$, since we need a dust correction prior to doing our inference. Dust corrections are particularly non trivial for stars with poor parallax signal to noise; these stars have a likelihood that is consistent with infinite distance and can therefore have severe dust corrections. This is most obvious in the giant stars, which tend to have the lowest signal to noise, and can therefore have dust corrections $>2$ mag. Although the dust model is probabilistic, as well as our distance inference, to move forward, we break from a fully Bayesian framework in which we might infer the distance and the dust simultaneously (more about this in the discussion). Instead, we correct for dust using a point estimate from the 3D dust map (Green et al., 2015). We determine this dust correction by iteratively inferring, and then sampling, our posterior PDF over $\varpi_{\mathrm{true}}$. We first generate our prior using the observed, attenuated photometry. Using this raw prior, we infer more precise parallaxes to all the stars in our sample. We use this more precise parallax posterior to get a measurement of dust in the 3D dust map for each star. We then iterate this process ten times, however the dust values seem to converge after a few iterations. The iterative process allows for slight variations in distance and corresponding dust to become more consistent with the CMD.

In detail, we take the $5\%$ quantile of each parallax posterior (the closest part) and query the 3D dust map at that distance. We take the $5\%$ quantile—as apposed to the median—to do a minimal dust correction because over-correcting for dust can push the likelihood function into an unphysical, extremely blue part of the CMD and thus bias the parallax inference. At the 3D position of each star, we sample the mode=median value of the probabilistic dust map. We apply each dust correction to our 2MASS photometry using Equation 12, where $Q_{\lambda}=[0.709,0.302]$ for bands $[J,K_{s}]$ respectively (see Table 6 in Schlafly & Finkbeiner 2011). The dust values we obtain in our final iteration, and which we apply to the photometry in our inference, are visualized in galactic coordinates in Figure 3. After applying the dust correction, we do not update the uncertainties, which is technically wrong and under-represents the true uncertainties; this leads to a less severe deconvolution, as mentioned previously. In other words, the features of our inferred CMD are conservatively broad.333This can be understood by considering the following procedure: if one incrementally increased the noise levels (leaving the data unchanged; just increasing the declared uncertainties), the deconvolved distribution would be increasingly narrower, since it targets the distribution that would generate the noiseless observations, as well as the observed data once noise is added. Thus, under-estimating the uncertainties leads to a conservative under-deconvolution. We could add the covariance of the dust to our likelihoods, but this is beyond the scope of this demonstration of concept.

4.2 The empirical prior

Using XD and the MMLE method described in the methods section, we generate the prior, shown in Figure 4. The left panel shows 1.4 million samples from the prior distribution, and the right panel shows $1\sigma$ contours for each component in the underlying Gaussian mixture model. Compared to the raw data in Figure 2, the deconvolved color–magnitude diagram for the Gaia + 2MASS stars is tighter, as expected. We’ll come back to more features in the data below, especially those seen in the posterior distribution.

4.3 Shrinkage of parallax uncertainties

Figure 5 shows the various probability distributions of our inference for a few randomly-chosen objects: the likelihood (black), prior (green), and posterior PDF (blue) on the true parallaxes. For the majority of the Gaia stars, our prior slightly increases the precision of the posterior PDF compared with the Gaia likelihood. Our parallax posterior PDFs are more precise than the TGAS catalog entries by a median factor of $1.2$. Examples of these slight increases in precision are visualized in Examples 1, 3, and 4 in Figure 5 which show slightly narrower posterior PDFs compared to the likelihoods. However, $14\%$ of our parallax posterior PDFs are more precise than the TGAS catalog by a factor of $>2$. This larger change in precision is visualized in Example 5, which shows a significantly narrower posterior PDF compared to the likelihood. The nature of our prior is also visible: it is apparent that it is made of a mixture of components, and that its impact, and shape, strongly depend on the location of each star in CMD space.

Figure 6 shows the distribution of the expectation values of the posterior parallax PDFs projected back onto the CMD space. The left panel shows the distribution for all 1.4 million stars in our dataset, and the right panel shows a subsample of the dataset and includes the $1\sigma$ uncertainties on the posterior PDFs. Compared with the raw data shown in Figure 2, it is clear the precision of the posterior PDFs is greater, especially for the red giant branch stars.

This change in precision is further illustrated in Figure 7, where we show the fractional changes in the variance. The left panel shows the natural log of the fractional change in variance as a function of $(J-K_{s})^{C}$ color, with the $1\sigma$ and $2\sigma$ contours of the distribution over plotted. The points are colored by the natural log of their fractional change in variance (the $y$ axis of the left panel) in both panels to help guide your eye for the right panel. The regions of CMD space with the greatest improvements are shown in black, and the regions with the lowest improvements are shown in yellow. The giants show obvious large improvements, as well as some main sequence stars. Regions of large improvement are those with the largest photometric and parallax errors, which are the faintest observed objects. This is a fairly standard result which we illustrated in the toy model above: the deconvolution of uncertainties is stronger in noisier regions of the data, leading to narrower features in the upper part of the CMD, for example. Stars with a color $(J-K_{s})^{C}\sim 0.6$ show a decrease in precision due to the width of the prior there; the more vertical structure in the red giant branch of the CMD.

Figure 8 shows the cumulative distribution function of the natural log of the fractional change in variance of the parallax posterior pdf relative to the TGAS catalog (the y-axis of Figure 7). This more directly shows what was alluded to in Figure 5. For the majority of the Gaia stars, our prior slightly increases the precision of the posterior PDF compared with the TGAS variance. Our parallax posterior PDFs are more precise than the TGAS catalog entries by a median factor of $1.2$, shown as the dotted line at $y=0.5$. The precision of our posterior PDFs is greater than or equal to the precision of the TGAS catalogue for $83\%$ of the stars, shown as the dotted line at $x=0$, and $14\%$ of our parallax posterior PDFs are more precise than the TGAS catalog by a factor of $>2$, which corresponds to a fraction decrease in variance by a factor of $>4$, shown as the dotted line at $y=0.14$.

We have shown that the posterior parallax estimates are more precise than the likelihoods, both on average and especially for low signal-to-noise stars. It is also interesting to ask whether our data-driven prior leads to biased inferences. In Figure 9, we show the shifts of the parallax estimates from likelihood TGAS catalog entry to posterior expectation. This shows that while many parallax estimates change substantially, they do not change on average, even at low signal-to-noise ratios. Our conclusion is that the prior developed here does not introduce additional bias into the TGAS parallax measurements.

4.4 Data products

The method produces a mixture-of-Gaussian parallax posterior PDF for every star in the intersection of TGAS, 2MASS, and the dust map. Associated with this paper is an electronic table containing our results444http://voms.simonsfoundation.org:50013/8kM7XXPCJleK2M02B9E7YIYmvu5l2rh/ServedFiles/, and everything that is necessary to use them responsibly. It includes the TGAS identifier, sky position, and dust-corrected photometry. It contains posterior expectations and variances for the stellar parallax. It also contains the full posterior PDF for each star, in the form of a set of PDF values, evaluated on a common grid in parallax. Because all of the data are public, and all of the code used for this project is available publicly under an open-source license, power users can generate any of the figures in this paper or any other data outputs straightforwardly.

4.5 Distances to M67

As a test of the accuracy of these photometric parallax inferences, we inferred the parallaxes to stars in the open cluster M67. M67 is at a Heliocentric distance of $800$–$900~{}{\rm pc}$ with an estimated age of $3$–$5~{}{\rm Gyr}$. Previous work (Bailer-Jones, 2015) has noted that distance estimates directly derived from the inverse of the observed (noisy) parallax are unreliable (significantly biased and uncertain) when the parallax signal-to-noise (SN) is below 5, $\varpi/\sigma_{\varpi}<5$.

With TGAS parallaxes having a minimum uncertainty of about $0.3~{}{\rm mas}$, the SN reaches 5 for the most certain parallax measurements at about $700~{}{\rm pc}$, just before the distance to this cluster. M67 is therefore a great test case for comparing parallax inferences with or without a prior.

We select sources from our catalog within a $1~{}{\rm deg}^{2}$ window centered on the position of M67; the cluster is visible as an over-density of stars in the sky positions alone. To be clear, we do no prior selection on previously known members of M67, but include all stars observed by TGAS within this window on the sky. Keep in mind, the (low) number of M67 members seen here is constrained by the (unknown) selection effects of the TGAS dataset. With our converged dust values, and its associated prior, we calculate the posterior over $\varpi_{\mathrm{true}}$ for each star within this window on the sky around M67. Figure 10 shows the comparison of the likelihoods and posteriors. The vertical red lines bracket the previously measured parallaxes of, or distances to, the cluster. The likelihoods of the observed parallaxes, shown in the top row in black, are broad and the cluster is not at all obvious. The posterior PDFs, shown in the bottom row in blue, are more sharply peaked showing an increase in precision. The posterior PDFs also have modes that are more similar to the previously determined distance to the cluster (vertical lines) for far more stars. This shows that our prior is not only making parallaxes more precise but also (possibly) increasing the accuracy.

5 Discussion

We have shown that it is possible to obtain photometric parallaxes for distant stars in the Gaia TGAS catalog without any use of physical stellar models, nor stellar density models of the Milky Way. We used the geometric parallaxes to calibrate a photometric model that is purely statistical, which is a model of the data rather than a model of stars per se. This opens up the possibility of completing the goals of the Gaia Mission without building in unnecessary assumptions about the mechanical properties of stars or the Galaxy.

We obtained the photometric parallaxes in this project by building a data-driven model of CMD of the stars in the TGAS data set, and using it as a prior PDF for Bayesian inference. The posterior PDFs for distance that we obtain are, in general, much narrower than the likelihood functions delivered by the Gaia Mission, and therefore the distance estimates (or, equivalently, parallax estimates) are much more precise. It is not surprising that a Bayesian inference provides more precise inferences than the likelihood function alone; Bayesian inferences bring in new information that decrease variance (but can introduce bias). We have shown that, in addition to pure precision improvements, at least some aspects of accuracy have been improved as well: The posterior distance estimates to stars in (at least one) stellar cluster are much more clustered than likelihood-based distance estimates, and they cluster around a sensible value. Accuracy is demonstrated by this test because each star in the cluster is subject to a different, unique distance prior (see Figure 5) and yet the posterior PDFs are all pulled to the same parallax value.

Although we have not performed principled hierarchical Bayesian inference in this project, we built our prior for each star’s parallax by performing a statistically responsible deconvolution of the Gaia TGAS data that is justifiable under a clear set of assumptions. This deconvolution, using a reasonable noise model and an assumption of stationarity, shrinks the parallax uncertainties, most dramatically for the stars measured at lower signal-to-noise. Because the method we use, XD, which is an Empirical Bayesian maximum-marginalized-likelihood estimator, accounts for heteroskedastic noise, we are able to build our photometric parallax model—which is a model of the CMD—using all the data, not just the data with the highest signal-to-noise. So, the generated model for the CMD is representative for our subset of the TGAS Catalog. Because our model is representitive of our full subset catalog, we expect any (sensible) parallax estimates we generate from our parallax posterior PDFs to be only weakly biased.

5.1 Comparison with other priors

There are other options for prior information to include when inferring distances to stars. One option is the exponentially declining spatial density (EDSD) of stars (Astraatmadja & Bailer-Jones, 2016a). It is parameterized by a scale length, and has the nice property that it is smooth out to very large distances so that the posteriors are also very smooth. It is a fairly weak prior. Figure 11 shows a comparison of the EDSD prior with our CMD prior. Here the parallax posterior PDF expectation value, and $1\sigma$ uncertainty, for a subset of stars in our dataset is converted into an absolute magnitude and visualized as the CMD. The left panel visualizes results using the EDSD prior, and the right panel visualizes results using our CMD prior. The variance for the CMD prior is significantly smaller than the variance for the EDSD prior. In general, EDSD is a fairly broad prior so the median Gaia star has a posterior variance equal to the likelihood variance. EDSD has the largest effect on the most noisy and negative parallaxes, and has the nice property of bringing these noisy or negative measurements to reasonable parallaxes.

5.2 What’s That Feature?

The denoised expectation values of the parallax projected into the CMD show many interesting and familiar features; we visualize this CMD again in Figure 12 and highlight some of these features.

The lower main sequence is very narrow while the upper main sequence is much wider. This may be an age and metallicity variation on the main sequence mixed with a Malmquist bias of seeing a much larger volume of upper main sequence stars (Malmquist, 1922). There is a plausible binary sequence of stars with about twice the brightness as the main sequence. The helium burning red clump is prominent for the volume DR1 probed, and the red giant branch is fairly straight and narrow. The main sequence turn off is noticeable and also surprisingly narrow, possibly a reflection of the star formation history of the local volume also mixed with a Malmquist bias.

5.3 Critical discussion of assumptions

In Section 3 we listed a set of assumptions upon which the photometric parallax method in this project is based. Here we return to these assumptions, and discuss them critically in more detail, under the same set of assumption labels:

stationarity

We have assumed that (locally) the morphology of the CMD is independent of signal-to-noise; that is, we have assumed that both the low and high signal-to-noise objects are being drawn from the same distribution. This isn’t precisely true because nearby stars in TGAS are mostly main sequence stars, while we can see luminous giants to much larger distances. So the red giant branch is a population of relatively lower signal-to-noise objects compared with the main sequence. This does not violate the stationarity assumption in detail, because the low and high signal-to-noise objects have to be drawn from the same distribution within some feature or patch of the CMD space. Thus, this has a small effect on our results, and could be improved by modeling the noise distribution simultaneously with the CMD.

selection

The biggest limitation of our method for generating the CMD from the data is that we have not used a TGAS selection function, completeness estimate, or inverse-selection-volume corrections. This model is a model of the TGAS Catalog (or really the TGAS–2MASS–PS1 intersection), not of the stars in the Milky Way, nor any volume-limited subsample of the Milky Way. Importantly, because we used all of the stars in the TGAS–2MASS–PS1 intersection, the model is representative of the full intersection, even the distant parts, where all the parallaxes are measured at low signal-to-noise. This is in contrast to techniques that might build the model from only the high signal-to-noise sources; such a model would be biased towards stellar types and compositions found locally, as well as lower-luminosity stars, and stars which happened to get good parallaxes. However, our use of the entire TGAS–2MASS–PS1 intersection without any selection function restricts the use of our CMD model as a prior to inferences within that same TGAS–2MASS–PS1 intersection. It would be biasing to use this same prior to infer distances for the full billion-star catalog released in Gaia DR1 alongside TGAS, or for stars fainter or brighter or bluer than those included in the TGAS–2MASS–PS1 intersection. A similar argument can be made for internal selection effects, such as inhomogeneities of the depth and noise distributions, which we have neglected. Those inhomogeneities are unlikely to affect the CMD results; we effectively model the average distribution of the TGAS–2MASS–PS1 intersection. By modeling the selection function in detail, one could improve the quality of the inferred CMD, its effect on individual objects, and its connection to physical models of the Milky Way.

big data

The empirical Bayesian methodology is a good approximation to full hierarchical inference in the limit of large numbers of data points, such that no individual star is carrying a lot of weight in the model of the CMD. With its millions of stars, the TGAS Catalog appears to safely live in this limit. However, because there are parts of the CMD that are poorly populated, even in TGAS, the empirical Bayes approximation may be bad for some portions of the CMD. In particular, there appear to be essentially no white dwarf stars, and other classes (like bright red giant branch stars) have few to no exemplars in the data set at high signal-to-noise (in parallax). For all these reasons, we do think that this approximation is not perfectly safe, and it is a goal to explore full hierarchical modeling in the coming years (and see, for example, Leistedt & Hogg 2017).

noise model

Everything in both the establishment of the data-driven prior and the use of it for parallax inference has been done under the deep assumption that the Gaia noise model is correct, and an additional assumption that the parallax uncertainty is always significantly larger than the photometric uncertainty in the $J$ band, even after dust correction. In the unlikely event (given these choices) that the uncertainties are over-estimated, we could be over-deconvolving. That is, the data-driven model could be producing a tighter or more informative distribution than what would truly be observed in a much higher signal-to-noise experiment or survey. That is, incorrect uncertainty estimates will bias the empirical-Bayes prior we obtain. On the inference side, wrong uncertainty estimates lead to further biases in the derived photometric parallaxes. Since we are explicitly ignoring the (non-Gaussian, in parallax space) photometric and dust-extinction uncertainties, our noise estimates are certainly biased, and this will translate into small biases in posterior parallax estimates.

dust

Within our probabilistic framework, we improperly handle dust. We are taking a point estimate of the probabilistic dust model at a quantile of our distance posterior, correcting our photometry using this point estimate, and then re-deriving the prior. A full, proper account would sample the dust model at distances sampled from the posterior, and propagate those samples back into the prior during the iterative process. The reason we don’t do a full probabilistic treatment, is that this proper dust distribution would not be a Gaussian distribution, which XD requires. We could calculate this more complex dust distribution and estimate it as a Gaussian distribution to feed to XD, however, this is beyond the scope of this paper.

mixture of Gaussians

The CMD model is a mixture of Gaussians, and furthermore fixed at $K=128$ components. This setting was chosen after some experimentation as being able to capture features but still easily optimized. We performed some experiments that suggest that larger $K$ would lead to better models, even in the conservative sense of a cross-validation score. Choosing an optimized or optimal number for this mixture is a natural extension of this work.

Beyond this, the choice of the Gaussian form itself for the mixture is questionable. The Gaussian has great properties, that we capitalize on both at the prior-generation stage (with XD), and at the photometric-parallax stage (capitalizing on Gaussian product rules). We could move to other forms, but they would be extremely expensive, computationally.

no physics

Although it is presented as an advantage of our method that we never make any use of physical models, it is also brings some significant disadvantages: For example, the CMD contains features that are obvious products of stellar evolution, but doesn’t in itself inform the theory or models of stellar evolution, because our mixture-of-Gaussian CMD model is divorced from any theoretical ideas about stars. For another, stellar evolution models are in fact very predictive and successful; our model does not capitalize in any way on those successes. A next-generation model might be designed to be data-driven but nonetheless capitalize on physics-driven successes. One framework for that would be to fit not the CMD but rather deviations of the CMD away from a good physics-based model prediction. This would also require modeling the selection function of the data under consideration.

The outputs of this project include posterior PDFs for stellar parallaxes. It is tempting to treat these outputs as equivalent to some kind of catalog of measurements, as we treated the original parallax measurements from TGAS. That is, it is tempting to treat these as simply “better” measurements. Taken one star at a time, this is permitted and true, and the basis for our claim that we have de-noised the TGAS Catalog. However, there are important differences, as there are with all probabilistic catalogs (Hogg & Lang 2011; Portillo et al. 2017), between the original TGAS measurements and the de-noised measurements. For users of our output, the most important difference between our output and the TGAS data input is that the TGAS data is the representation of a likelihood function, whereas our ouptut is a representation of a posterior PDF, one per star. Posterior information used naively can lead to serious statistical errors, because each data point has had a prior multiplied in. Multiplicative uses of the data will effectively take the prior PDF to a large power. We have (in previous work) given examples of correct uses of posterior samplings or PDFs for subsequent inferences (Hogg et al. 2010; Foreman-Mackey et al. 2014); we encourage power users to consult those methodological contributions before engaging with these outputs.

6 Conclusion

We forgo an explicit conclusion; we attempt to summarize the conclusions in the Abstract above.

It is a pleasure to thank Ana Bonaca (Harvard), Andy Casey (Monash), Stephen Feeney (Flatiron), Dustin Lang (Toronto), David Spergel (Flatiron), and the attendees at the Stars Group Meeting at the Flatiron Institute Center for Computational Astrophysics for comments and input. This project was developed in part at the 2016 NYC Gaia Sprint, hosted by the Center for Computational Astrophysics at the Simons Foundation in New York City. This work has made use of data from the European Space Agency (ESA) mission Gaia555http://www.cosmos.esa.int/gaia, processed by the Gaia Data Processing and Analysis Consortium666http://www.cosmos.esa.int/web/gaia/dpac/consortium (DPAC). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research was partially supported by the NSF (AST-1517237), NASA (grant NNX12AI50G), and the Moore-Sloan Data Science Environment at NYU. BL was supported by the National Aeronautics and Space Administration through Einstein Postdoctoral Fellowship Award Number PF6-170154. It made use of the NASA Astrophysics Data System. All the code used in this project is available online777https://github.com/andersdot/photoParallax under an open-source license.