Hierarchical probabilistic inference of the color-magnitude diagram and shrinkage of stellar distance uncertainties

TL;DR

This paper introduces a hierarchical probabilistic model that enhances stellar distance estimates by integrating noisy parallax and photometric data through a data-driven color-magnitude diagram, improving accuracy especially for uncertain measurements.

Contribution

It presents a novel hierarchical model that infers a noise-deconvolved color-magnitude diagram without relying on stellar models, significantly improving distance estimates for Gaia stars.

Findings

Distance SNR improved by 19% on average

Uncertainty reduced by 36% on average

8% of objects have SNR doubled

Abstract

We present a hierarchical probabilistic model for improving geometric stellar distance estimates using color--magnitude information. This is achieved with a data driven model of the color--magnitude diagram, not relying on stellar models but instead on the relative abundances of stars in color--magnitude cells, which are inferred from very noisy magnitudes and parallaxes. While the resulting noise-deconvolved color--magnitude diagram can be useful for a range of applications, we focus on deriving improved stellar distance estimates relying on both parallax and photometric information. We demonstrate the efficiency of this approach on the 1.4 million stars of the Gaia TGAS sample that also have APASS magnitudes. Our hierarchical model has 4~million parameters in total, most of which are marginalized out numerically or analytically. We find that distance estimates are significantly improved for the noisiest parallaxes and densest regions of the color--magnitude diagram. In particular, the average distance signal-to-noise ratio and uncertainty improve by 19~percent and 36~percent, respectively, with 8~percent of the objects improving in SNR by a factor greater than 2. This computationally efficient approach fully accounts for both parallax and photometric noise, and is a first step towards a full hierarchical probabilistic model of the Gaia data.

Tables

Table 1. Table 1: Summary of our notation.

$s$	object index (the $s$-th star)
$d_s,{\varpi }_s,M_s,C_s$	true distance, parallax, absolute magnitude, and color
${\hat{\varpi }}_s,{\sigma }_{{\hat{\varpi }}_s}^2$	parallax estimate and its variance
${\hat{m}}_s,{\hat{C}}_s,{\sigma }_{{\hat{m}}_s}^2,{\sigma }_{{\hat{C}}_s}^2$	apparent magnitude and color estimates, and their variances
$b_s$	index of the color–magnitude bin of the $s$th object
$b$	generic index of color–magnitude bin (the $b$th bin)
$n_b$	object count in the $b$-th color–magnitude bin
$\left\{n_b\right\}$	set of all galaxy counts $n_b$, summing to $N_{\mathrm{stars}}$
$f_b$	fractional galaxy count in the $b$-th color–magnitude bin
$\left\{f_b\right\}$	set of all fractional bin counts $f_b$, summing to $1$
$\{d_s,b_s\}$	distances and bins of all stars in the sample
$\{{\hat{m}}_s,{\hat{C}}_s\}$	all magnitude and color estimates

Equations

$$\displaystyle p(d_{s}|\hat{m}_{s},\hat{C}_{s},\hat{\varpi}_{s})=\int{\mathrm{d}}M_{s}\ {\mathrm{d}}C_{s}\ p\bigl{(}\hat{m}_{s},\hat{C}_{s},\hat{\varpi}_{s}\bigr{\rvert}M_{s},d_{s},C_{s}\bigr{)}\ p\bigl{(}M_{s},d_{s},C_{s}\bigr{)}.$$

$$p\left(\hat{\varpi}_{s}\bigr{\rvert}d_{s}\right)=\mathcal{N}\bigl{(}\hat{\varpi}_{s}-d_{s}^{-1};\sigma_{\hat{\varpi}_{s}}^{2}\bigr{)},$$

$$p\bigl{(}\hat{m}_{s},\hat{C}_{s}\bigr{\rvert}M_{s},d_{s},C_{s}\bigr{)}=\mathcal{N}\bigl{(}M_{s}+5\log_{10}d_{s}-\hat{m}_{s};\sigma_{\hat{m}_{s}}^{2}\bigr{)}\ \mathcal{N}\bigl{(}\hat{C}_{s}-C_{s};\sigma_{\hat{C}_{s}}^{2}\bigr{)}.$$

$$p\left(M_{s},C_{s}\bigr{\rvert}\{f_{b}\}\right)=\sum_{b=1}^{B}f_{b}\ K_{b}(M_{s},C_{s}),$$

$$p(d_s$$

$$\displaystyle\quad\quad\times\ \$$

$$\displaystyle p\left(b_{s}\bigr{\rvert}\bigl{\{}f_{b}\bigr{\}}\right)\$$

$$\displaystyle p\left(M_{s},C_{s}\bigr{\rvert}b_{s}\right)\$$

$$\displaystyle p\left(\bigl{\{}f_{b}\bigr{\}}\bigr{\rvert}\bigl{\{}d_{s},b_{s},\hat{m}_{s},\hat{C}_{s},\hat{\varpi}_{s}\bigr{\}}\right)\ =\ p\bigl{(}\bigl{\{}f_{b}\bigr{\}}\bigr{\rvert}\{n_{b}\}\bigr{)}\ \propto\ \prod_{b}\frac{f_{b}^{n_{b}}}{n_{b}!},$$

$$\displaystyle p\left(b_{s}\bigr{\rvert}\bigl{\{}f_{b}\bigr{\}},\hat{m}_{s},\hat{C}_{s},\hat{\varpi}_{s}\right)\$$

$$g_{b_s}$$

Full text

stars: distances, stars: C–M diagrams, methods: statistical.

\journalinfo

Prepared for ApJ

1 Introduction

Gaia and LSST will soon provide parallaxes for many millions of stars inhabiting our Galaxy (Gaia Collaboration et al., 2016; LSST Science Collaboration et al., 2009). However, parallax measurements for most of these objects are very noisy due to parallax errors significantly increasing with distance. Most analyses adopt signal-to-noise cuts, thus discarding large amounts of hard-won data. In rare cases, uncertainties are explicitly accounted for, via probabilistic hierarchical models, for example (e.g.,  Sesar et al., 2016). Regardless of how parallax or distance uncertainties are accounted or ignored, reducing them can improve a variety of studies and increase the scientific return of a mission like Gaia.

In addition to direct geometric determination via parallaxes, stellar distances can be indirectly inferred from spectral information. This is because the vast majority of stars live in compact regions of color–magnitude space, such as the main sequence or the red giant branch. This information can be used as a prior for disentangling the contributions of distance and absolute magnitude in apparent magnitudes. In other words, photometric information can complement parallax information and prove crucial to determine distances. The influence of distance priors, as well as color–magnitude priors based on stellar models, was explored in previous work (e.g.,  Astraatmadja & Bailer-Jones, 2016b, a), which relied on theoretical models or simulations. Yet, stellar models predict narrow tracks in color–magnitude, which are reductive of the diversity of real objects, and may introduce biases when used as priors for determining distances. In this paper, we address this problem and present a purely data-driven approach for improving parallax-based distance estimates with color–magnitude information inferred directly from the data, with no use of external models.

Our approach relies on a standard property of hierarchical probabilistic models: by using a large number of noisy observations (here, parallaxes and magnitudes) one can construct an estimate of the underlying distribution (here, a color–magnitude diagram), which is in turn applied as a prior to shrink the uncertainties on the properties of the individual objects (here, the distances). This Bayesian shrinkage naturally occurs within hierarchical models when computing or sampling the full joint posterior distribution. In other words, the whole is greater than the sum of the parts, and stellar distance estimates can be improved by simultaneously analyzing all the available data. This is particularly relevant for Gaia, which will soon provide many millions of sources with parallax, proper motion, and photometric measurements, offering the opportunity to construct accurate models of the distribution of magnitudes, colors, dust, and stellar densities in the Galaxy, with minimal external data. Those models will in turn significantly improve our estimates of the intrinsic properties of individual stars. The core goal of this paper is to demonstrate how one aspect of this hierarchical inference can be done in practice. We focus on evaluating the shrinkage of distance estimates obtained by constructing a flexible color–magnitude diagram based on the joint analysis of millions of noisy parallaxes and apparent magnitudes. We make a number of simplifying assumptions to make our demonstration more transparent, such as adopting uninformative distance priors and point estimates of dust corrections. Those assumptions could be relaxed or complexified in subsequent work.

In parallel to this work, a noise-deconvolved color–magnitude diagram reconstruction from the Gaia data with a variant of the extreme deconvolution (XD) algorithm Bovy et al. (2011) is presented in Anderson et al. (2017). Five key elements distinguishes the XD inference from the present work. First, the XD model is based on 2MASS photometry, which is more precise than the APASS data we use here. Second, in the XD model, the color–magnitude diagram is modeled as a (transformed) Gaussian mixture where the positions and sizes of the mixture components are also inferred. Third, in the XD model, the noise model is more approximate (magnitude uncertainties are ignored) to make the method more tractable. Fourth, in the XD model dust is a smaller correction (because of the infrared photometry) and it is iteratively optimized to be self-consistent with the final color–magnitude diagram model. Fifth, the XD model is a maximum-marginalized likelihood method; it is not fully probabilistic at all levels of the hierarchy in the hierarchical model. Those differences allow to construct a color–magnitude diagram at higher resolution, at the cost of an approximate methodology scheme (Anderson et al., 2017). The methodology and results presented here fully account for a realistic noise model (for both parallaxes and photometry), but they rely on noisier data and keep the resolution (the positions and widths of the Gaussian mixture component) fixed, resulting in a more reliable inference, but of a less precise color–magnitude diagram. Suitably combining the two approaches will be the key to fully exploiting the full-Mission Gaia data with computationally tractable but nonetheless justified hierarchical probabilistic models.

This paper is structured as follows: our model and inference framework are described in Sec. 2, and applied to the Gaia data in Sec. 3. Conclusions and perspectives are presented in Sec. 4.

2 Model

We consider a set of stars indexed as $s=1,\cdots,{N_{\rm stars}}$, each characterized by a (latent) distance $d_{s}$, absolute magnitude $M_{s}$, and color $C_{s}$. We consider only one magnitude and one color for simplicity, but it should be noted that the model and method presented below can be straightforwardly extended to multiple magnitudes and colors.

Intrinsic properties like distance and absolute magnitude are not observable, they must be inferred from apparent magnitudes and parallax measurements. The estimate of the parallax is denoted $\hat{\varpi}_{s}$ and is assumed to have a Gaussian variance $\sigma_{\hat{\varpi}_{s}}^{2}$. We will consider two magnitudes, $\hat{m}_{s}$ and $\hat{m}^{\prime}_{s}$, assumed to be uncorrelated and to have Gaussian variances $\sigma_{\hat{m}_{s}}^{2}$ and $\sigma_{\hat{m}^{\prime}_{s}}^{2}$. We will use the first one $\hat{m}_{s}$ as a reference magnitude for inferring the absolute magnitude $M_{s}$, and the second one to form a color estimate $\hat{C}_{s}=\hat{m}^{\prime}_{s}-\hat{m}_{s}$ with Gaussian variance $\sigma_{\hat{C}_{s}}^{2}=\sigma_{\hat{m}_{s}}^{2}+\sigma_{\hat{m}^{\prime}_{s}}^{2}$.

We aim at estimating the distance $d_{s}$ of each star from the noisy data $\hat{m}_{s}$, $\hat{C}_{s}$ and $\hat{\varpi}_{s}$. While distance is directly connected to the parallax via $\varpi_{s}=1/d_{s}$, it is also informed by the apparent magnitude since $m_{s}=M_{s}+5\log_{10}d_{s}$ where $d_{s}$ is expressed in units of $10$ pc. Note that when only the apparent magnitude is available, distance and absolute magnitude are degenerate and cannot be disentangled. This degeneracy is partially broken with the parallax information. Here, we seek to incorporate the knowledge that stars do not have arbitrary colors and magnitude. The way this information enters distance estimates is made obvious by writing the posterior probability distribution on the distance,

[TABLE]

This integral marginalizes over the true absolute magnitude and color. This might be expensive to perform numerically, but the choices we will make below will allow us to execute it analytically.

The first term of Eq. (1) is a likelihood function, and the second term is the prior. Assuming that the magnitude and parallax estimates are independent, the likelihood function factorizes as the product of two terms,

[TABLE]

and

[TABLE]

The final term, $p(M_{s},d_{s},C_{s})$, is the prior knowledge about the distances, magnitudes, and colors of stars. It is typically factorized as $p(M_{s},C_{s}|d_{s})p(d_{s})$, where $p(d_{s})$ is the distance prior, which could depend on position and based on a three-dimensional model of the Galaxy. The color–magnitude prior is typically made distance-independent, $p(M_{s},C_{s})$, and based on stellar models (e.g.,  Astraatmadja & Bailer-Jones, 2016b, a).

In this work, we will adopt a uniform distance prior and focus on the magnitude–color term, which we parametrize as p(M_{s},C_{s}\bigr{\rvert}\{f_{b}\}). We construct a model of the relative abundance of objects in color–magnitude cells (i.e., in two dimensions). The color–magnitude distribution is described as a linear mixture of $B$ components,

[TABLE]

with $K_{b}$ the $b$th kernel. The parameters $\{f_{b}\}$ refer to the relative probabilities of finding objects in the various kernels, and must sum to one ($\sum_{b}f_{b}=1$).

For the kernels, we adopt Gaussian distributions to make the integral of Eq. (1) analytically tractable. The $b$-th kernel will be centered at $(\mu_{b,0},\mu_{b,1})$ and have a diagonal covariance $(\sigma_{b,0}^{2},\sigma_{b,1}^{2})$. We take $\mu_{b+1,0}-\mu_{b,0}=\sigma_{b,0}$ and $\sigma_{b,0}$ constant (similarly for the color dimension) to uniformly and contiguously tile a rectangular region of interest in color–magnitude space. With this parameterization, the integral of Eq. (1) is tractable and leads to

[TABLE]

Finally, to facilitate parameter inference, we will introduce a latent variable $b_{s}$ denoting the bin the $s$th object belongs to. Then, we can equivalently write the color–magnitude model as

[TABLE]

Our notation is summarized in Table 1. Note that in all of the above we have assumed that all magnitudes are properly dereddenned, i.e., that the absorption by interstellar dust has been corrected for. As discussed below, we will use the distance point estimate $1/\hat{\varpi}_{s}$ for performing this correction. While this is formally incorrect (since the reddenning correction depends on distance which is a parameter we are inferring), this approximation will be sufficient for the case-study of this paper.

2.1 Inference

Since we fix the kernel locations $\{(\mu_{b,0},\mu_{b,1})\}$ and covariances $\{(\sigma_{b,0}^{2},\sigma_{b,1}^{2})\}$, our color–magnitude model is fully described by the relative amplitudes $\{f_{b}\}$. Those could be set by prior knowledge (e.g., external data or stellar models), in which case one could use Eq. (2) to infer the distance of each object. Here, we seek to infer $\{f_{b}\}$ from the data. Thus, the full posterior of interest is $p(\{d_{s}\},\{f_{b}\}|\{\hat{m}_{s},\hat{C}_{s},\hat{\varpi}_{s}\})$, which has $B+{N_{\rm stars}}$ parameters.

Given the number of parameters and the natural degeneracies between magnitudes and distances, standard sampling techniques may be difficult to apply. However, the inference is greatly simplified by focusing on $p(\{b_{s}\},\{f_{b}\}|\{\hat{m}_{s},\hat{C}_{s},\hat{\varpi}_{s}\})$, where distances are marginalized over. We adopt a Gibbs sampling strategy (see e.g.,  Casella & George, 1992; Wandelt et al., 2004; Levin et al., 2009; Brooks et al., 2011), which consists of alternating draws from the conditional distributions of $\{b_{s}\}$ and $\{f_{b}\}$. At the $i$th iteration, we will draw new values of the $\{f_{b}\}$ and $\{b_{s}\}$ parameters given the values of the previous iteration. In other words, we first draw $\{f_{b}\}^{(i)}$ given $\{b_{s}\}^{(i-1)}$ (and the data), and for each object then draw $b_{s}^{(i)}$ given $\{f_{b}\}^{(i)}$ (and the data). The sequence $\{f_{b}\}^{(i)},\{b_{s}\}^{(i)}$ for $i=1,\cdots,N_{\mathrm{samples}}$ forms a Markov Chain with the target posterior distribution of interest as equilibrium distribution. This allows us to avoid the magnitude–distance degeneracies and exploit other properties of our model, such as conjugate priors. We now detail how to draw from the correct conditional distributions.

The first draw is fairly standard: with the bin locations $\{b_{s}\}$ fixed, the fractional amplitudes $\{f_{b}\}$ follow a Dirichlet distribution entirely determined by $\{n_{b}\}$, with $n_{b}$ the number of objects in the $b$-th bin. All the other parameters enter the constant proportionality factor, so the target conditional distribution is

[TABLE]

where we have assumed the simplest uninformative prior for the relative amplitudes, i.e., that they are positive and normalized to one, p(\bigl{\{}f_{b}\bigr{\}})=\delta^{D}(1-\sum_{b}f_{b})\prod_{b}\Theta(f_{b}), where $\delta^{D}$ is the Dirac delta function and $\Theta$ the step function. Given the number counts $\{n_{b}\}$, one can draw $\{f_{b}\}$ from this distribution using standard tools for Dirichlet distributions. Alternatively, one could draw $\gamma_{b}\sim\mathrm{Gamma}(n_{b}+1)$ and take $f_{b}=\gamma_{b}/\sum_{b}^{\prime}\gamma_{b}^{\prime}$.

The second step of the Gibbs sampler is to draw the bins $\{b_{s}\}$ given the amplitudes $\{f_{b}\}$. Those draws can be performed independently (thus, in parallel) over objects. This may not look simple since we have to marginalize over the distances, and the target conditional distribution to sample is

[TABLE]

However, closer examination reveals that the bins can be simply drawn from a multinomial distribution with amplitudes given by

[TABLE]

with

[TABLE]

Since we keep the kernel locations and covariances fixed, those weights can be precalculated via ${N_{\rm stars}}$ independent dimensional integrals. We evaluate them numerically given that the integrand is simple and the integration can be tuned easily. If this operation were prohibitive, or if the kernel parameters were not kept fixed, one could sample the distance $d_{s}$ (jointly or conditionally on the bin $b_{s}$). Since the kernels and likelihoods admit simple gradients and Hessians, one could use Hamiltonian Monte Carlo (e.g.,  Duane et al., 1987; Neal, 2012) to efficiently sample all ${N_{\rm stars}}$ distances simultaneously.

2.2 Discussion

We now briefly discuss some of the advantages and limitations of our approach.

First, we restricted our attention to the color–magnitude diagram, and neglected any dependency on other quantities such as galactic latitude. We also adopted uniform distance priors. This can be corrected by adopting or constraining models of the 3D distribution stellar densities, described by kernel mixtures, for example (see e.g.,  Astraatmadja & Bailer-Jones, 2016b, a). Although this extension is technically trivial, we have not developed it since we focus on how color–magnitude information informs distance estimates without the use of stellar models. Similarly, our framework could be extended to other observables such as velocities.

Second, our kernel mixture model offers a significant amount of freedom for describing the color–magnitude diagram, but could be improved in various ways. Changing the kernels would not strongly affect our inference framework, unless they were not differentiable, in which case it would not possible to use Hamiltonian Monte Carlo. Analytic marginalization of true magnitudes and colors was made possible by adopting Gaussian kernels. Furthermore, we did not optimized the positions and covariances of the kernels, unlike in standard Gaussian Mixture models (e.g.,  Bovy et al., 2011; Anderson et al., 2017). Compared to those, our tiling of color–magnitude space requires more components (many of which are zero) but is easy to initialize, and also converges quickly.

Third, we assumed that the magnitudes are perfectly dust-corrected. However, dust extinction depends on distance, which is a parameter of our model. Furthermore, reliable 3D dust maps are only available for a limited region of three-dimensional space. In principle, dust corrections should really be inferred jointly with the absolute magnitudes and colors of the data at hand. This is explored in Anderson et al. (2017).

3 Application to Gaia TGAS

We consider the Gaia data (Gaia Collaboration et al., 2016), specifically the first data release (DR1) of the Tycho-Gaia astrometric solution (hereafter TGAS, Lindegren et al., 2016). We restrict our attention to the objects with valid B and V magnitudes from the AAVSO Photometric All Sky Survey (APASS) Data Release 9 (Munari et al., 2014; Henden & Munari, 2014). We also remove objects with parallax signal-to-noise ratio (SNR) lower than 1. This leads to 1.4 million objects with magnitude, color, and parallax information. We do not apply more stringent parallax or color cuts since the purpose of our method is exactly to construct a color–magnitude model from both low- and high-SNR objects. Finally, we apply dust corrections based on position and distance point estimate ($1/\hat{\varpi}$) with the three-dimensional dust map of Green et al. (2015). Our data sample is summarized in Fig. 1, which shows the magnitude, color, and parallax SNR distributions. The bulk of the objects has parallax SNR lower than 10 and is at $M_{V}<4$, in the upper part of the color–magnitude diagram. This highlights the need for a correct inference framework that exploits all objects, since focusing on high-SNR objects would bias the results and prevent us from correctly inferring the fainter regions of the color–magnitude space. This is illustrated below. Note that we did not add the recommended parallax error offset of 0.3 (Lindegren et al., 2016) because overestimated errors typically cause excessive deconvolution. For the type of hierarchical inference we are considering, underestimated errors lead to more conservative results.

We create a validation sample by extracting 10% of the objects, adding significant amount of noise to the parallax estimates, and verifying that our framework improves the distances consistently with the original values. Further details about this process are given below. We also split the main sample according to parallax SNR, into two samples of equal size containing the ‘best’ and ‘worst’ parallaxes. We perform the inference on those two samples as well as their union.

For each of the three samples, we use the Gibbs sampler presented above to draw thirty thousand samples of the fractional amplitudes and bins. Thus, our model has $B+{N_{\rm stars}}$ sampled parameters, and $3{N_{\rm stars}}$ parameters marginalized over (the true magnitudes, colors, and distances), either numerically or analytically. The positions and widths of the kernels are obtained by uniformly tiling (with 60 points) the color and magnitude lines in the ranges $[-0.3,1.7]$ and $[-1,7]$, respectively.

The top panels of Fig. 2 show the input data. In particular, the right panel shows a random sample of two thousand objects with their associated errors. Note that the absolute magnitude Gaussian error shown for each object is approximate and obtained by linearizing the parallax likelihood and adding the resulting error to the photometric one in quadrature. This is formally incorrect (and only adopted for visualization purposes) since it underestimates the true absolute magnitude error. Our framework does not make this approximation since it explicitly models those components.

The mean and standard deviation of the resulting color–magnitude diagram (with bins marginalized over) are shown in the bottom panels of Fig. 2. Those were obtained by computing the mean and standard deviation of the set of models obtained by applying Eq. (4) to all the posterior samples we drew. As expected, the models recovered are significantly narrower than the data since we are effectively deconvolving observational errors to produce a noiseless color–magnitude diagram. The main sequence, its turn-over, and the giant branch, are successfully recovered at higher resolution than in the input noisy data. Features corresponding to binary stars and red clump stars are not clearly visible due to the coarse binning we adopted. While we could increase the resolution of our model, there is no guarantee that those features will be resolved since the photometric noise is significant in our data. Adopting Gaia and 2MASS photometry could help resolve this issue (see Anderson et al., 2017). In addition, our Gibbs sampling scheme is sensitive to the resolution adopted, and the mixing of parameters slows down as the resolution increases. One could circumvent this obstacle by adopting a Hamiltonian Monte Carlo sampling scheme. This is possible since the gradients with respect to the parameters are simple and can be efficiently calculated. Finally, we note that the oscillatory features are also the result of an excessive resolution of the model in the noisiest regions of color–magnitude space. This could be avoided by adopting variable kernel size (inferred or tuned a priori).

Fig. 3 shows a comparison of the stellar distances with and without the extra information brought by the color–magnitude model. For each object we compute the distance posterior distribution via Eq. (1) on a grid, and we average over samples of the bin allocations and relative amplitudes of the color–magnitude model, which is equivalent to marginalizing them out. Our point of comparison is the set of parallax-only distance estimates, i.e., not using our hierarchical model but geometric information only, via Eq. (2). Note that neither the full nor the marginalized posterior distributions are Gaussian, as expected and also shown below. Nevertheless, the mean and standard deviation provides a useful metrics. Our inference methodology provides samples of the full, joint non-Gaussian posterior distribution, and those samples should be used for any subsequent inference relying on the color–magnitude diagram to correctly capture the data and model uncertainties. We measure that on average the distance SNR improves by 19% and the distance uncertainty decreases by 36%. We find that 8% of the objects have their distance uncertainty halved after the inclusion of color–magnitude information via our hierarchical model. This shrinkage of the distance uncertainties is most efficient in the densely populated regions of color–magnitude space. Fig. 4 shows the distance posterior distributions obtained with our method for a few randomly chosen objects, illustrating the uncertainty shrinkage. We also place these objects in the color–magnitude diagram; the errors on the absolute magnitude are obtained by resampling. Note that as in Fig. 2, this underestimates the absolute magnitude uncertainty, but nevertheless provides a qualitatively useful error estimate.

Fig. 5 shows the mean and standard deviation of the color–magnitude diagram resulting from performing the inference on the subsamples with parallax SNR cuts (i.e., splitting our main sample at parallax SNR of 5). Those demonstrate that including the noisiest objects is essential for correctly inferring the faint regions of magnitude space. The main sequence is well recovered with the high-SNR objects, while the red giant branch is barely detected. By contrast, it is well recovered with the low-SNR objects, but the main sequence is then partially erased. This is a natural consequence of the SNR increasing with absolute magnitude. This highlights the importance of a correct probabilistic framework, capable of correctly exploiting data with heterogeneous noise to reconstruct the noiseless color–magnitude diagram. To illustrate the effect of the SNR cut we also resampled the data and kept objects that satisfied the cuts. This provides a first-order estimate of the detection probability or selection function of those two subsamples, shown in the right panels of Fig. 5, illustrating what regions of the color–magnitude diagram are unexplored.

We now turn to the validation sample. Since we do not know the true distances for most objects in our data set, we add noise to the parallax estimate, at a level equal to ten times the parallax error. We then compute the posterior distribution on the distance (on a distance grid) with and without color–magnitude information, as in Fig. 4. We simply use the mean model shown in Fig. 2 as a color-magnitude prior. The results are shown in Fig. 6. Given that those objects have significant amounts of noise, causing the distance posterior distribution to be highly non-Gaussian, the mean distance overestimates the true distance (the original parallax-based estimate). The hierarchical model significantly decreases this effect, i.e., improves the distance estimates both in terms of mean and uncertainty, demonstrating the validity of our inference scheme.

Finally, we perform an additional test of our method on open clusters. We retrieve the coordinates and distances of known open clusters from the WEBDA database111www.univie.ac.at/webda. From our TGAS-APASS sample, we select the objects within a radius corresponding to 3 pc around each cluster, using the known cluster distance. We only keep objects with parallaxes consistent with the open cluster distance at the three sigma level, to reduce the number of candidate members. Fig. 6 shows the result of applying our framework to those objects for a few selected clusters, i.e., using the color–diagram inferred above to inform the distance estimates. Even though not all objects selected by this procedure are true cluster members, most distance estimates undergo a visible improvement toward the cluster distance, with the SNRs mildly improved. In addition, Fig. 4 showed that the distance posterior distributions produced by the hierarchical model are significantly more Gaussian, i.e., the typical heavy tail produced by the parallax-only likelihood function for low SNR objects is strongly decreased. This approach could potentially help identifying cluster members and determining cluster distances with minimal external information.

4 Conclusions

Stellar distances are ubiquitous in astronomy and are most directly determined geometrically via parallaxes. However, parallax measurements are unreliable for the faintest, most distant objects, which inhabit most of the Galaxy. We presented a framework for obtaining improved distance estimates with both parallax and color–magnitude information without external data or priors such as stellar models. This exploits a well-known property of probabilistic hierarchical models: by utilizing the wealth of information concealed in millions of data points (the noisy and less noisy ones), one can model the underlying distributions and in turn improve the quality of the individual parameter estimates. We presented a version of this idea where the color–magnitude distribution of stars is directly inferred from noisy parallax and magnitude measurements via an efficient parameter inference scheme. This is the first color–magnitude model consistently capturing both parallax and photometric uncertainties. We applied this methodology to Gaia+APASS data and demonstrated that it leads to significant improvements in the distance estimates, specifically for the faintest, most distant objects.

The framework described here included a number of simplifying assumptions. First, we used uniform distance priors and described the color–magnitude diagram as a mixture of Gaussian kernels with fixed positions and widths. Those two assumptions could easily be relaxed without affecting our inference methodology or the quantitative results presented here. Second, we ignored selection effects that may distort the color–magnitude diagram and the distance distribution. Those could be incorporated but require external data and modifications to the sampling method. Third, we neglected the dependency of the color–magnitude on other properties (e.g., galactic latitude), and we performed the magnitude reddenning corrections using existing three-dimensional dust maps evaluated at the parallax-based distance point estimates. Our framework could also support those extensions with moderate technical changes. However, in order to fully take advantage of the available data, one should jointly infer the three-dimensional distribution of dust and stellar density, and compare with results based on physical models and priors. Any selection effects modulating the detection probability or the measured properties or uncertainties of stars should also be included. This requires more substantial technical developments due to the large number of interlaced parameters. We intend to explore this avenue in future work.

We thank Lauren Anderson, Andy Casey, Jo Bovy, Adrian Price-Whelan, Semyeong Ho, Keith Hawkins, David Spergel, and the Stars Group Meeting at the Flatiron Institute, for stimulating discussions and useful comments on this manuscript. BL was supported by NASA through the Einstein Postdoctoral Fellowship (award number PF6-170154). DWH was partially supported by the NSF (AST-1517237) and the Moore–Sloan Data Science Environment at NYU. 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 Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund.