seestar: Selection functions for spectroscopic surveys of the Milky Way

TL;DR

This paper introduces a Bayesian Poisson likelihood method for calculating and combining selection functions in spectroscopic surveys of the Milky Way, enabling better understanding of observational biases and intrinsic stellar properties.

Contribution

It presents a generalizable, robust approach for determining and transforming selection functions applicable to various multi-object spectrographs and survey overlaps.

Findings

Successfully recreates synthetic spectroscopic samples from mock catalogues.

Handles overlapping survey fields and combined samples effectively.

Transforms sky-based selection functions into intrinsic stellar property functions.

Abstract

Selection functions are vital for understanding the observational biases of spectroscopic surveys. With the wide variety of multi-object spectrographs currently in operation and becoming available soon, we require easily generalisable methods for determining the selection functions of these surveys. Previous work, however, has largely been focused on generating individual, tailored selection functions for every data release of each survey. Moreover, no methods for combining these selection functions to be used for joint catalogues have been developed. We have developed a Poisson likelihood estimation method for calculating selection functions in a Bayesian framework, which can be generalised to any multi-object spectrograph. We include a robust treatment of overlapping fields within a survey as well as selection functions for combined samples with overlapping footprints. We also provide a method for transforming the selection function that depends on the sky positions, colour, and apparent magnitude of a star to one that depends on the galactic location, metallicity, mass, and age of a star. This `intrinsic' selection function is invaluable for chemodynamical models of the Milky Way. We demonstrate that our method is successful at recreating synthetic spectroscopic samples selected from a mock galaxy catalogue.

Tables

Table 1. Table 1: Parameters and hyperparameters of the GMMs and their respective priors to be fit to the photometric and spectroscopic catalogues. i 𝑖 i refers to the field on which the model fit is performed with every field fit independently.

Data	Component	Model	Parameters	Prior	Hyperparameters
Photometric	Distribution Function	GMM	$N_i$	$U[0,inf]$
			${\pi }_{ik}$	Dirichlet	$\alpha =1/K$
			$w_{ik}=N_i{\pi }_{ik}$
			${𝝁}_{ik}$, ${𝚺}_{ik}$	NIW	${𝐦}_0=\mathrm{med}[{𝐯}_{\min },{𝐯}_{\max }]$, $\lambda ={10}^{-4}$,
					$𝚿={\left(\frac{\mathrm{diag}\left({𝐯}_{\max }-{𝐯}_{\min }\right)/2}{5}\right)}^2$, $\nu =2$
Spectroscopic	Distribution Function	GMM	$N_i$	Poisson	$\mathrm{mean}=N_i^{\mathrm{phot}}$
			${\pi }_{ik}$	Dirichlet	${\alpha }_{ik}^{\prime }$
			$w_{ik}=N_i{\pi }_{ik}$
			${𝝁}_{ik}$, ${𝚺}_{ik}$	NIW	${𝐦}_{0ik}^{\prime }$, ${\lambda }_{ik}^{\prime }$
					${𝚿}_{ik}^{\prime }$, ${\nu }_{ik}^{\prime }$
	Selection Function	GMM	$\widetilde{{\kappa }_{ik}}$	$U[-inf,inf]$
			$\widetilde{w_{ik}}=\sqrt{\left|2\pi \widetilde{{𝚺}_{ik}}\right|}{\mathrm{logit}}^{-1}\left(\widetilde{{\kappa }_{ik}}\right)$
			$\widetilde{{𝝁}_{ik}}$, $\widetilde{{𝚺}_{ik}}$	NIW	$\widetilde{{𝐦}_0}=\mathrm{med}[{𝐯}_{\min },{𝐯}_{\max }]$, $\tilde{\lambda }={10}^{-4}$,
					$\widetilde{𝚿}={\left(\frac{\mathrm{diag}\left({𝐯}_{\max }-{𝐯}_{\min }\right)/2}{5}\right)}^2$, $\tilde{\nu }=2$

Table 2. Table 2: Notation followed in this paper.

$𝐱=(𝜽,𝐯)$	Full coordinates of the star.
$\left\{X_{\mathrm{phot}}\right\}\left(\left\{X_{\mathrm{spec}}\right\}\right)$	Sample of photometric (spectroscopic) data.
$𝜽=(l,b)$	Galactic coordinates.
$𝐯=(m,c)$	Observable coordiantes: colour and apparent magnitude.
$\tau ,\left[M/H\right],{ℳ}_{\mathrm{ini}},s$	Intrinsic coordinates: age, metallicity, initial mass and distance
${\widetilde{ℳ}}_{\mathrm{ini}}$	Initial mass scaled to the range ${\widetilde{ℳ}}_{\mathrm{ini}}\in [0,1]$.
${\Theta }_i$	Event that a star is on field $i$.
$S$	Event that a star is selected in the spectroscopic catalogue.
$S_i$	Event that a star is selected on field $i$ of the spectrosopic catalogue.
$f_{\mathrm{phot}}\left(f_{\mathrm{spec}}\right)$	Distribution function (spectroscopic distribution).
$n_{\mathrm{phot}}\left(n_{\mathrm{spec}}\right)$	Density of stars in the photometric (spectrosopic) sample.
$N_{\mathrm{phot}}\left(N_{\mathrm{spec}}\right)$	Number of stars in the photometric (spectrosopic) sample.
$g_{\mathrm{DF}}\left(g_{\mathrm{SF}}\right)$	GMM for the distribution function (selection function)
$ϵ\left(\widetilde{ϵ}\right)$	Parameters of the distribution function (selection function).
${\alpha }_k$	Dirichlet concentration parameters for the photometric density GMM.
$w_{ik},{𝝁}_{ik},{𝚺}_{ik}$	Parameters of field $i$, component $k$ of the photometric density GMM.
${\pi }_{ik}$	Normalised component loadings of the photometric density GMM.
$N$	Normalisation of the photometric density GMM.
$\widetilde{w_{ik}},\widetilde{{𝝁}_{ik}},\widetilde{{𝚺}_{ik}}$	Parameters of field $i$, component $k$ of the selection function GMM.
$\widetilde{{\kappa }_{ik}}$	logit transform of the component weights of the selection function GMM.
${ℒ}^{\mathrm{phot}}\left({ℒ}^{\mathrm{spec}}\right)$	Likelihood function for the photometric density (spectroscopic density) fit.
$M$	Number of fields or patches.
$K$	Number of GMM components in photometric density.
$\tilde{K}$	Number of GMM components in selection function.

Equations

$$\mathbb{P}(\mathrm{S}\mid \mathbf{x})=\frac{\mathbb{P}(\mathbf{x}\mid \mathrm{S})\mathbb{P}(\mathrm{S})}{\mathbb{P}(\mathbf{x})}.$$

$$\mathbb{P}(\mathrm{S}\mid \mathbf{x})=\frac{f^{\mathrm{spec}}(\mathbf{x})\mathbb{E}\left[N^{\mathrm{spec}}\right]}{f^{\mathrm{phot}}(\mathbf{x})\mathbb{E}\left[N^{\mathrm{phot}}\right]}=\frac{\mathbb{E}\left[n^{\mathrm{spec}}(\mathbf{x})\right]}{\mathbb{E}\left[n^{\mathrm{phot}}(\mathbf{x})\right]},$$

$$\mathbb{P}({\Theta }_i\mid \mathbf{\theta })=\{\begin{array}{ll}1& \text{if θ is located on patch i},\\ 0& \mathrm{otherwise}.\end{array}$$

$$\mathbb{P}({\mathrm{S}}_i\mid \mathbf{\theta },\mathbf{v})=\mathbb{P}({\mathrm{S}}_i\mid {\Theta }_i,\mathbf{v})\mathbb{P}({\Theta }_i\mid \mathbf{\theta })$$

$$\mathbb{P}(\mathrm{S}\mid\boldsymbol{\theta},\mathbf{v})=\mathbb{P}\left(\bigcup_{i=1}^{M}\mathrm{S}_{i}\,\bigg{|}\,\boldsymbol{\theta},\mathbf{v}\right)$$

$$\begin{split}\mathbb{P}\left(\bigcup_{i}^{M}\mathrm{S}_{i}\,\bigg{|}\,\boldsymbol{\theta},\mathbf{v}\right)=\sum_{k=1}^{M}(-1)^{k+1}\left[\sum_{1\leq i_{1}<...<i_{k}\leq M}\mathbb{P}(\mathrm{S}_{i_{1}},\mathrm{S}_{i_{2}}...\mathrm{S}{i_{k}})\right]\end{split}$$

$$\begin{array}{r}\mathbb{P}({\mathrm{S}}_1,{\mathrm{S}}_2,...)=\prod _{i=1,2,...}\mathbb{P}({\mathrm{S}}_i\mid {\Theta }_i,\mathbf{v})\mathbb{P}({\Theta }_i\mid \mathbf{\theta }).\end{array}$$

$$\mathbb{P}({\mathrm{S}}_i\mid \mathbf{v},{\Theta }_i)=\frac{\mathbb{P}(\mathbf{v}\mid {\mathrm{S}}_i,{\Theta }_i)\mathbb{P}({\mathrm{S}}_i\mid {\Theta }_i)}{\mathbb{P}(\mathbf{v}\mid {\Theta }_i)}.$$

$$\mathbb{P}({\mathrm{S}}_i\mid {\Theta }_i,\mathbf{v})=\frac{f_i^{\mathrm{spec}}(\mathbf{v}\mid {\Theta }_i)\mathbb{E}\left[N_i^{\mathrm{spec}}\right]}{f_i^{\mathrm{phot}}(\mathbf{v}\mid {\Theta }_i)\mathbb{E}\left[N_i^{\mathrm{phot}}\right]}=\frac{\mathbb{E}\left[n_i^{\mathrm{spec}}(\mathbf{v}\mid {\Theta }_i)\right]}{\mathbb{E}\left[n_i^{\mathrm{phot}}(\mathbf{v}\mid {\Theta }_i)\right]},$$

$$G(\mathbf{v}\mid \mathbf{\mu },\mathbf{\Sigma })=\frac{1}{\sqrt{\mid 2\pi \mathbf{\Sigma }\mid }}\exp \left(-\frac{(\mathbf{v}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{v}-\mathbf{\mu })}{2}\right),$$

$$n_i^{\mathrm{phot}}(\mathbf{v}\mid {\mathbf{ϵ}}_i,{\Theta }_i)=\sum _{k=1}^K{w}_{ik}G(\mathbf{v}\mid {\mathbf{\mu }}_{ik},{\mathbf{\Sigma }}_{ik}),$$

$$\begin{array}{rl}\ln ({\mathcal{L}}^{\mathrm{phot}}(X_i^{\mathrm{phot}}|{\mathbf{ϵ}}_i))=-\int & \mathrm{d}\mathbf{v}{n}_i^{\mathrm{phot}}(\mathbf{v}\mid {\mathbf{ϵ}}_i,{\Theta }_i)\\ & +\sum _{j=1}^{N_i^{\mathrm{phot}}}\log \left(n_i^{\mathrm{phot}}({\mathbf{v}}_j\mid {\mathbf{ϵ}}_i,{\Theta }_i)\right),\end{array}$$

$$\begin{array}{rl}\ln ({\mathcal{L}}^{\mathrm{phot}}(X_i^{\mathrm{phot}}|{\mathbf{ϵ}}_i))\propto -& N+N_i^{\mathrm{phot}}\ln (N)\\ & +\sum _{j=1}^{N_i^{\mathrm{phot}}}\ln \left({\hat{n}}_i^{\mathrm{phot}}({\mathbf{v}}_j\mid {\mathbf{ϵ}}_i,{\Theta }_i)\right),\end{array}$$

$$\begin{array}{rl}\mathrm{NIW}(& \mathbf{\mu },\mathbf{\Sigma }\mid {\mathbf{m}}_0,\lambda ,\mathbf{\Psi },\nu )\\ \propto & \mathcal{N}\left(\mathbf{\mu }\mid {\mathbf{m}}_0,\frac{1}{\lambda }\mathbf{\Sigma }\right)|\mathbf{\Sigma }{|}^{-\frac{(\nu +p+1)}{2}}\exp \left[-\frac{1}{2}\mathrm{Tr}\left(\mathbf{\Psi }{\mathbf{\Sigma }}^{-1}\right)\right],\end{array}$$

$$\mathrm{DIR}({\pi }_1,{\pi }_2...{\pi }_K|{\alpha }_1,{\alpha }_2...{\alpha }_K)\propto \prod _{k=1}^K{\pi }_k^{{\alpha }_k-1}$$

$$\mathrm{BIC}=\ln (N_i^{\mathrm{phot}})d-2\ln ({\mathcal{L}}_{\max }^{\mathrm{phot}})$$

$$n_i^{\mathrm{spec}}(\mathbf{v}\mid {\Theta }_i)=n_i^{\mathrm{phot}}(\mathbf{v}\mid {\Theta }_i)\times \mathbb{P}({\mathrm{S}}_i\mid \mathbf{v},{\Theta }_i).$$

$$g_i^{\mathrm{SF}}(\mathbf{v}\mid \widetilde{{\mathbf{ϵ}}_i},{\Theta }_i)=\sum _{k=1}^{\tilde{K}}\widetilde{w_{ik}}G(\mathbf{v}\mid \widetilde{{\mathbf{\mu }}_{ik}},\widetilde{{\mathbf{\Sigma }}_{ik}}),$$

$$\begin{array}{rl}\ln ({\mathcal{L}}^{\mathrm{spec}}(X_i^{\mathrm{spec}}|& {\mathbf{ϵ}}_i,\widetilde{{\mathbf{ϵ}}_i}))\propto \\ -\int \mathrm{d}\mathbf{v}& n_i^{\mathrm{phot}}(\mathbf{v}\mid {\mathbf{ϵ}}_i,{\Theta }_i)g_i^{\mathrm{SF}}(\mathbf{v}\mid \widetilde{{\mathbf{ϵ}}_i},{\Theta }_i)\\ +& \sum _{j=1}^{N_i^{\mathrm{spec}}}\log \left(n_i^{\mathrm{phot}}({\mathbf{v}}_j\mid {\mathbf{ϵ}}_i,{\Theta }_i)g_i^{\mathrm{SF}}({\mathbf{v}}_j\mid \widetilde{{\mathbf{ϵ}}_i},{\Theta }_i)\right).\end{array}$$

$$\begin{array}{rl}\ln (\mathbb{P}({\mathbf{ϵ}}_i,\widetilde{{\mathbf{ϵ}}_i}|X_i^{\mathrm{spec}}))=\ln (& {\mathcal{L}}^{\mathrm{spec}}(X_i^{\mathrm{spec}}|{\mathbf{ϵ}}_i,\widetilde{{\mathbf{ϵ}}_i}))\\ -& N+N_i^{\mathrm{phot}}\ln (N)\\ +& \mathrm{DIR}({\pi }_1,{\pi }_2...{\pi }_K|{\alpha }_1^{\prime },{\alpha }_2^{\prime }...{\alpha }_K^{\prime })\\ +& \sum _{k=1}^K\left[\mathrm{NIW}({\mathbf{\mu }}_{ik},{\mathbf{\Sigma }}_{ik}|{\mathbf{m}}_{0ik}^{\prime },{\lambda }_{ik}^{\prime },{\mathbf{\Psi }}_{ik}^{\prime },{\nu }_{ik}^{\prime })\right]\\ +& \sum _{k=1}^{\tilde{K}}\left[\mathrm{NIW}(\widetilde{{\mathbf{\mu }}_{ik}},\widetilde{{\mathbf{\Sigma }}_{ik}}|\widetilde{{\mathbf{m}}_0},\tilde{\lambda },\tilde{\mathbf{\Psi }},\tilde{\nu })\right].\end{array}$$

$$\mathrm{BIC}=\ln (N_i^{\mathrm{spec}})d-2\ln ({\mathcal{L}}_{\max }^{\mathrm{spec}})$$

$$\mathbb{P}(\mathrm{S}\mid s,[M/H],{\mathcal{M}}_{\mathrm{ini}},\tau ,\mathbf{\theta }),$$

$$\begin{array}{r}\mathbb{P}(\mathrm{S}\mid s,\tau ,[M/H],{\mathcal{M}}_{\mathrm{ini}},\mathbf{\theta })=\mathbb{P}(\mathrm{S}\mid \mathbf{v}(s,\tau ,[M/H],{\mathcal{M}}_{\mathrm{ini}}),\mathbf{\theta }).\end{array}$$

$$m=M+5{\log }_{10}\left(\frac{s}{10\mathrm{pc}}\right)$$

$$(c,M)=F^{\mathrm{iso}}(\tau ,[M/H],{\mathcal{M}}_{\mathrm{ini}}),$$

$${\tilde{\mathcal{M}}}_{\mathrm{ini}}=\frac{{\mathcal{M}}_{\mathrm{ini}}-{\mathcal{M}}_{\mathrm{ini}}^{\min }}{{\mathcal{M}}_{\mathrm{ini}}^{\max }-{\mathcal{M}}_{\mathrm{ini}}^{\min }},$$

$$l,b=\{\begin{array}{l}30,60(\text{Field 1}),\\ 30,61(\text{Field 2}).\end{array}$$

$$\begin{array}{rl}\mathbb{P}(S|m_H,J-K)=0.1& \times \left(1-\tanh \left(\frac{m_H-13.5}{e^{-2}}\right)\right)/2\\ & \times \left(1-\tanh \left(\frac{J-K-0.5}{e^{-2}}\right)\right)/2.\end{array}$$

$$m_x=M_x+5{\log }_{10}\left(\frac{s}{10\mathrm{pc}}\right)+A_x(l,b,s),$$

$$\begin{array}{r}\mathbb{P}(\mathrm{S}\mid \mathbf{\theta },\mathbf{v})=\mathbb{P}({\mathrm{S}}_{\mathrm{A}}\mid \mathbf{\theta },\mathbf{v})+\mathbb{P}({\mathrm{S}}_{\mathrm{B}}\mid \mathbf{\theta },\mathbf{v})-\mathbb{P}({\mathrm{S}}_{\mathrm{A}},{\mathrm{S}}_{\mathrm{B}}\mid \mathbf{\theta },\mathbf{v}).\end{array}$$

Full text

We have developed a Poisson likelihood estimation method for calculating selection functions in a Bayesian framework, which can be generalised to any multi-object spectrograph. We include a robust treatment of overlapping fields within a survey as well as selection functions for combined samples with overlapping footprints. We also provide a method for transforming the selection function that depends on the sky positions, colour, and apparent magnitude of a star to one that depends on the galactic location, metallicity, mass, and age of a star. This ‘intrinsic’ selection function is invaluable for chemodynamical models of the Milky Way. We demonstrate that our method is successful at recreating synthetic spectroscopic samples selected from a mock galaxy catalogue.

††pubyear: 2018††pagerange: seestar: Selection functions for spectroscopic surveys of the Milky Way–12

1 Introduction

Within our Milky Way, we are able to resolve stars out to large distances, but with that privilege comes a bias. The most comprehensive spectroscopic surveys can typically measure quantities for only a fraction of a percent of the stars in the Milky Way. The survey or ‘observed’ selection function is the probability of a star being included in the survey given its sky positions, colour and apparent magnitude. To make inferences regarding the intrinsic chemodynamical structure of the Milky Way, knowledge is required of the intrinsic selection function, i.e. the probability of a star being included in the survey given its intrinsic coordinates: Galactic location, metallicity, mass, and age.

The second data release (Gaia Collaboration et al., 2018) of the ESA’S Gaia mission has pushed us far further than ever before, measuring photometry, positions, and proper motions for over 1.3bn stars. In parallel, several ground-based spectroscopic surveys are measuring spectra for millions of these stars. Despite an often relatively simple nominal selection in colour and apparent magnitude, taking cross-matches with other surveys or selecting stellar type sub-samples result in a selection function that is no longer simple. Furthermore, to convert the observed colour-magnitude selection function into an intrinsic one that depends on distance, metallicity, mass, and age, one needs to engage with stellar isochrones. In order to understand the bias generated by selecting subsamples, we must calculate the selection functions for these surveys.

Many studies have developed survey selection functions, in which the completeness along a line of sight is given by the ratio of the number of stars in the spectroscopic survey to the number of stars in a photometric survey in the same region of colour and apparent magnitude. Das & Binney (2016) and Das et al. (2016) use this method to construct selection functions for halo blue horizontal branch stars and K giants in Sloan Extension for Galactic Understanding and Exploration-2 (SEGUE-II, Xue et al., 2011). Similar methods are used in determining the selection function for the Radial Velocity Experiment (RAVE, Kordopatis et al., 2013) by Wojno et al. (2017) and for the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST, Cui et al., 2012) Spectroscopic Survey of the Galactic Anticentre (LSS-GAC, Xiang et al., 2017) by Chen et al. (2018). Stonkutė et al. (2016) perform a similar analysis of the Gaia-ESO survey (Gilmore et al., 2012), taking into account the survey’s observing strategy.

Vickers & Smith (2018) generate a combined selection function for LAMOST (Zhao et al., 2012), RAVE (Kunder et al., 2017) and the Tycho-Gaia Astrometric Solution (TGAS, Michalik et al., 2015) by binning the fields in colour, apparent magnitude, and distance using a synthetic galaxy catalogue as the assumed complete sample. Bovy et al. (2012) and Bovy et al. (2014) consider the dependence of the selection function on apparent magnitude for G-type dwarfs in the Sloan Extension for Galactic Understanding and Exploration (SEGUE, Yanny et al., 2009) survey and the Apache Point Observatory Galaxy Evolution Experiment (APOGEE, Majewski et al., 2017). Their selection function along any line of sight is assumed to be uniform in apparent magnitude with limits defined either by the faintest star observed in a field or by the survey’s nominal magnitude limit. Bovy et al. (2014) also present a selection function in distance that takes into account the dust extinction along the line of sight.

Finally, Nandakumar et al. (2017) and Mints & Hekker (2019) determine selection biases for a large number of spectroscopic surveys by binning the sample in colour-magnitude space, either using a regular grid or a specialised median binning algorithm, and comparing with the 2-Micron All Sky Survey (2MASS, Skrutskie et al., 2006).

Amongst all previous works, fundamental aspects of the methodology remain the same. The dependence of the selection function on colour and magnitude is treated as a ratio of number counts of stars between the spectroscopic and a given photometric survey or a synthetically generated population, which is assumed to be complete. Uncertainties in measurements in colour and apparent magnitude are not used in the calculation. The effects of overlapping coordinate fields on the selection function are not considered. Selection functions are largely constructed as a function of colour and magnitude and not converted to intrinsic coordinates, with the exception of Das & Binney (2016), Das et al. (2016), and Sanders & Binney (2015). With the exception of Vickers & Smith (2018), no methods are presented for combining selection functions.

In this paper, we build on previous work to create seestar111https://github.com/aeverall/seestar, a Python code that can be applied to any multi-fibre spectroscopic survey, independent of its footprint on the sky, the number of stars observed, and the selection criteria of the survey. We construct an algorithm to treat the limitations of Poisson noise when bins have small numbers of stars by using a Poisson point process, whose parameters we determine using maximum likelihood estimation. We propose how the uncertainties in colour and magnitude measurements may be incorporated into this. We include a method for calculating the union of overlapping field probabilities, which also enables selection functions of independent surveys to be combined. We use isochrones to convert the selection function depending on colour and apparent magnitude to one that depends on distance, metallicity, mass, and age. This is an essential component of chemodynamical models of the Milky Way (Schönrich & Binney, 2009; Sanders & Binney, 2015; Das & Binney, 2016).

In Section 2, we demonstrate how to calculate the selection function in observable coordinates (colour and apparent magnitude) and intrinsic coordinates (distance, metallicity, mass, and age), and how to determine the dependence on sky position. The results of tests on a mock catalogue are presented in Section 3. Finally, we discuss our results and potential future developments and applications of this method in Section 4. Section 2 is rather technical in nature, and we therefore advise readers only interested in the performance of the selection function package to skip to the results.

Appendix C provides a table of parameter descriptions for reference along with a graphical Bayesian network to help with following our notation and understanding the method.

2 Method

The selection function $\mathbb{P}(\mathrm{S}\mid\mathbf{x})$ of a stellar survey is the probability of a star being in the survey, $\mathrm{S}$, given the star’s coordinates $\mathbf{x}$. The coordinates may be observed (sky positions, apparent magnitude, and colour) or intrinsic quantities (galactic location, metallicity, mass, and age). We assume that for the spectroscopic survey, there is a photometrically selected catalogue of stars which is complete in the region of colour and apparent magnitude observable by the spectrograph. Throughout this paper, we will use the superscript ‘spec’ to refer to the spectroscopic sample and ‘phot’ to refer to the photometric sample. By Bayes’ theorem

[TABLE]

$\mathbb{P}(\mathbf{x})\mathrm{d}^{n}x\equiv f(\mathbf{x})\mathrm{d}^{n}x$ is the probability that a star chosen at random has coordinates within the volume element, $\mathrm{d}^{n}x$, where $f(\mathbf{x})$ is the distribution function of stars in the Milky Way. $\mathbb{P}(\mathbf{x}\mid\mathrm{S})\mathrm{d}^{n}x\equiv f^{\mathrm{spec}}(\mathbf{x}\mid\mathrm{S})\mathrm{d}^{n}x$ is the probability that a star picked at random from the survey has coordinates in $\mathrm{d}^{n}x$, where $f^{\mathrm{spec}}(\mathbf{x})$ is the distribution function of observed stars. Finally, $\mathbb{P}(\mathrm{S})=\mathbb{E}\left[N^{\mathrm{spec}}\right]/\mathbb{E}\left[N\right]$ is the probability of a star entering the survey and is given by the ratio of the expected number of stars in the survey to the expected total number of stars in the Milky Way. Under the assumption that the photometric sample is complete within the given range of observable coordinates, we denote $f^{\mathrm{phot}}(\mathbf{x})=f(\mathbf{x})$ and $N^{\mathrm{phot}}=N$ and substitute into Equation (1).

[TABLE]

where $\mathbb{E}\left[n^{\mathrm{spec}}(\mathbf{x})\right]$ and $\mathbb{E}\left[n^{\mathrm{phot}}(\mathbf{x})\right]$ are the expected number densities of stars in the spectroscopic and photometric surveys respectively, at coordinates $\mathbf{x}$.

The observation of a star is dependent on the star’s observable properties, $\mathbf{x}=(l,b,c,m)$ where $(l,b)$ are the star’s coordinates on the sky and $(c,m)$ the colour and apparent magnitude of the star. For ease of notation we will group these into positional coordinates, $\boldsymbol{\theta}=(l,b)$ and photometric properties, $\mathbf{v}=(c,m)$.

The positional coordinates indicate which region of the sky or ‘patch’ the star belongs to. The best method to characterize the dependence of the selection function on these patches depends on the survey design. This is described in Section 2.1. On each patch, the selection function is calculated as a function of $\mathbf{v}$ as described in Section 2.2. In Section 2.3, the observed coordinates are transformed into intrinsic coordinates, $\mathbf{v}(s,\mathrm{[M/H]},\mathcal{M}_{\mathrm{ini}},\tau)$, where $s,\mathrm{[M/H]},\mathcal{M}_{\mathrm{ini}},\tau$ are distance, metallicity, initial mass, and age respectively. These coordinates underpin the description of the chemodynamical distribution of stars in the Milky Way.

2.1 Dependence of selection function on sky positions

To determine the dependence of the selection function on positional coordinates, we bin the sky into independent regions called patches across which the selection function does not vary directly as a function of sky position. Therefore we can denote $\Theta_{i}$ as the event that the star’s coordinates are located on patch $i$. Using this notation, we can define the probability of a star being located on a field given the star’s sky positions as

[TABLE]

2.1.1 Defining a patch

We assume that the selection function given $\mathbf{v}$ is independent of $\boldsymbol{\theta}$ across any patch. This motivates us to divide the sky into smaller patches particularly in locations where we expect rapid changes in the selection function with respect to sky positions. However, using smaller patches reduces the number of stars per patch, which amplifies the effects of Poisson noise. A sufficient coarseness is required such that the effects of Poisson noise are limited but fine enough such that the dependence of the selection function on sky positions is small.

Before continuing, we should briefly explain what we mean by Poisson noise and why it is so central to this work. Detections and recordings of stars in stellar catalogues can be considered as independent events which are randomly sampled from a distribution over observed coordinates. Thus the detection of stars by stellar surveys is well modelled by a Poisson point process. The signal-to-noise ratio for a Poisson point process scales with the square root of the mean signal strength, $\frac{\mathrm{S}}{\mathrm{N}}\sim\sqrt{\lambda}$ where $\lambda$ is the expected number of events per interval of observed coordinate space. Therefore the larger the number of events per interval, the stronger the signal-to-noise ratio. This is why reducing the number of stars per patch in order to increase the spatial resolution also leads to a reduction in the signal-to-noise ratio.

A multi-fibre spectroscopic survey is constructed by placing fibres on a plate so that each fibre is at the image of a star. The plate covers a solid angle on the sky. The solid angle of sky observed by each plate is here referred to as a field, and the coordinates of the centre of the field define the field pointing. The distributions of fields for APOGEE and RAVE in sky positions are shown in Figure 1 where each coloured circle represents the extent of a different field, and the field pointings are located at the centre of each circle. For multi-fibre spectroscopic surveys we define a patch as the region of sky covered by a field. From here on we will refer to patches as fields.

For a point on a single field, the selection function is the probability that the star will be observed by that field given the star’s photometric properties, $\mathbf{v}$. If a coordinate does not lie in a field, the selection function is [math].

The selection function for field $i$ is given by

[TABLE]

where $\mathrm{S}_{i}$ is the event that a star is selected by field $i$.

In many spectrograph surveys, as can be seen for APOGEE and RAVE in Figure 1, fields may heavily overlap leading to a single coordinate being observed multiple times on different fields. The selection function is the probability that at least one field selects the star. This is given by the union of the event of each field selecting the star.

[TABLE]

where $M$ is the total number of fields employed. This is calculated using all observations made on field $i$, even if the same star is observed on another field. The probability of the union of being on either of the two fields is the probability that one or the other occurs. We expand this in terms of the selection by individual fields

[TABLE]

where

[TABLE]

assuming the events $\mathrm{S}_{1},\mathrm{S}_{2}...$ are independent. Appendix A provides a detailed explanation of this expansion.

As mentioned earlier, many surveys contain multiple observations of the same positional coordinates. These can be observations taken on separate days to different magnitude depths. For Equation (7) to be appropriate, different observations of the same field should only be considered as separate fields if the observations are independent, i.e. if the probability of a star being selected by one observation is independent of whether the star is selected by the other observation. If the observations are dependent, as is the case if stars in one observation are chosen deliberately to be exclusive or inclusive of those observed in another observation, then the observations should be combined to form a single field in the selection function. In our method we combine all observations with the same field pointing as a single field.

In the following sections, we examine the dependence of the selection function on $\mathbf{v}$. These are calculated for a given field, $i$. Having calculated the selection function for each field as a function of $\mathbf{v}$, they are combined using the method above to achieve the full selection function for the entire survey.

2.2 Selection function in observable coordinates

The probability of a star being selected by field $i$ given $\mathbf{v}$, and given that the star lies on the field (i.e. $\mathbb{P}(\Theta_{i}\mid\boldsymbol{\theta})=1$) is

[TABLE]

$\mathbb{P}(\mathbf{v}\mid\Theta_{i})\mathrm{d}\mathbf{v}\equiv f(\mathbf{v}\mid\Theta_{i})\mathrm{d}\mathbf{v}$ is the probability that a star chosen at random on field $i$ has coordinates within the $\mathrm{d}\mathbf{v}$ volume element, where $f(\mathbf{v}\mid\Theta_{i})$ is the distribution of stars in the Milky Way inside the cone projecting onto field $i$.

$\mathbb{P}(\mathbf{v}\mid\mathrm{S}_{i},\Theta_{i})\mathrm{d}\mathbf{v}\equiv f^{\mathrm{spec}}(\mathbf{v}\mid\Theta_{i})\mathrm{d}\mathbf{v}$ is the probability that a star observed on field $i$ of the spectroscopic survey has coordinates within the $\mathrm{d}\mathbf{v}$ volume element. Finally, $\mathbb{P}(\mathrm{S}_{i})=\mathbb{E}\left[N_{i}^{\mathrm{spec}}\right]/\mathbb{E}\left[N_{i}\right]$ is the probability of a star on field $i$ entering the survey, and is given by the ratio of the number of stars in the survey on field $i$ to the total number of stars in the Milky Way inside the cone projecting onto field $i$. Assuming that the photometric sample is complete within the given range of observable coordinates $f^{\mathrm{phot}}(\mathbf{v}\,\mid\,\Theta_{i})=f(\mathbf{v}\,\mid\,\Theta_{i})$ and $N_{i}^{\mathrm{phot}}=N_{i}$. Substituting into Equation (8)

[TABLE]

where $n_{i}^{\mathrm{spec}}(\mathbf{v}\mid\Theta_{i})=f_{i}^{\mathrm{spec}}(\mathbf{v}\mid\Theta_{i})N_{i}^{\mathrm{spec}}$ is the number density of stars observed by the survey on field $i$ and $n_{i}^{\mathrm{phot}}(\mathbf{v}\mid\Theta_{i})=f_{i}^{\mathrm{phot}}(\mathbf{v}\mid\Theta_{i})N_{i}^{\mathrm{phot}}$ is the number density of stars in the Milky Way on the cone projecting onto field $i$.

2.2.1 Number density of photometric sample

We start by calculating the expected number density of stars in the Milky Way on field $i$, $\mathbb{E}\left[n_{i}^{\mathrm{phot}}(\mathbf{v}\mid\Theta_{i})\right]$. The choice of photometric survey is discussed in Section 4.

The stars in the photometric survey represent a Poisson realisation of the smooth underlying number density function, $n^{\mathrm{phot}}$. The aim is to use the observed stars to estimate the true smooth number density function. We assume this function can be parameterised as a bivariate Gaussian Mixture Model (GMM) for each field $i$. Each bivariate Gaussian component is given by

[TABLE]

where $\boldsymbol{\mu}\in\mathbb{R}^{2}$ is the mean of the bivariate Gaussian in colour and magnitude and $\boldsymbol{\Sigma}\in\mathbb{R}^{2\times 2}$ is the symmetric covariance matrix. Using this, we can write

[TABLE]

where $\boldsymbol{\epsilon}_{i}\in\mathbb{R}^{K\times 6}$ for the $K$ components of the GMM, each with six parameters defining the bivariate Gaussian component $w$, $\mu^{0}$, $\mu^{1}$, $\Sigma^{00}$, $\Sigma^{11}$ and $\Sigma^{01}$ ($\boldsymbol{\Sigma}$ is symmetric so $\Sigma^{10}=\Sigma^{01}$).

We then find the parameters $\boldsymbol{\epsilon}_{i}$ by maximising the likelihood of the photometric catalogue stars. The Poisson likelihood is derived from Poisson count probabilities to give

[TABLE]

where $\mathbf{v}_{j}$ are the colour and apparent magnitude of star $j$ in the photometric catalogue on field $i$ (see Appendix B).

We reparameterise the weights as normalised ‘loadings’, $\pi_{ik}=w_{ik}/N$ where $N=\sum_{k=1}^{K}\left[w_{ik}\right]$.

[TABLE]

where $\hat{n}^{\mathrm{phot}}_{i}=n^{\mathrm{phot}}_{i}/N$. The normalisation, $N$ is now a parameter of the model whilst the loadings provide $K-1$ free parameters under the constraint $\sum_{k=1}^{K}\pi_{ik}=1$.

To obtain the posterior probabilities on each of our parameters we need priors. The prior on the normalisation $N$ is uniform, $N\sim U\left[0,\infty\right]$. For the mean and covariance of the Gaussian components, we employ a Normal Inverse Wishart (NIW) prior.

[TABLE]

where $p$ is the number of degrees of freedom of the system, in our case two (colour and apparent magnitude). The Normal and Inverse Wishart distributions are the conjugate priors for the mean and covariance of a multivariate Gaussian distribution respectively. As a result, the posterior on $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ will also be represented by a NIW distribution. The chosen hyperparameters ($\mathbf{m}_{0},\lambda,\boldsymbol{\Psi},\nu$) are dependent on the prior information about the data. Choosing small values of $\lambda$ and setting $\nu=p$ provides a least informative prior on the parameters. We discuss a good choice of hyperparameters in our mock tests in Section 3.

The loadings are evaluated with a Dirichlet prior,

[TABLE]

which satisfies the constraint $\sum_{k=1}^{K}\left[\pi_{ik}\right]=1$. The Dirichlet distribution is the conjugate prior for membership probability of a categorical variable. Setting concentration hyperparameters $\alpha_{i}=\alpha=1/K\,\forall\,i$ provides a least informative prior.

For a given number of components, $K$ we determine the posteriors on parameters of the normalised GMM using the BayesianGaussianMixture module from scikit-learn (Pedregosa et al., 2011) with parameters initialised by $k$-means clustering. This module iterates using a Variational Inference algorithm which alternates between assigning component membership probabilities and fitting the individual components. We recover the best fit hyperparameters of the posterior Dirichlet and NIW distributions $\alpha_{k}^{\prime},\mathbf{m}_{0k}^{\prime},\lambda_{k}^{\prime},\boldsymbol{\Psi}_{k}^{\prime}$ and $\nu_{k}^{\prime}$ where $k$ is the Gaussian mixture component. These hyperparameters define the posterior distribution of parameters for the GMM. The posterior on the normalisation, $N$ is given by the first two terms of the right hand side of Equation (13) where the uniform prior doesn’t provide a contribution. This is proportional to a Poisson distribution with mean $N^{\mathrm{phot}}_{i}$.

To determine the number of components, $K$ we fit our GMM with $K=1$ up to a maximum of $K=20$ components and calculate the Bayesian information criterion (BIC, Schwarz, 1978) for the best fit parameters (taken as the expected values of the posteriors)

[TABLE]

where $d=K\times 6$ (the number of parameters in our model) and $\mathcal{L}^{\mathrm{phot}}_{\mathrm{max}}$ is the likelihood corresponding to the best fit parameters. Our final model for the photometric sample colour-apparent magnitude function is the model which minimizes the BIC.

2.2.2 Selection function in observed coordinates

The selection function in observed coordinates for field $i$ is $\mathbb{P}(\mathrm{S}_{i}\mid\mathbf{v},\Theta_{i})$. The stars in the spectroscopic catalogue represent a Poisson realisation of the product of the distribution function and the selection function

[TABLE]

Therefore, we can once more use the Poisson likelihood in Equation (12), replacing $n^{\mathrm{phot}}(\mathbf{v}\mid\Theta_{i})$ with $n^{\mathrm{spec}}(\mathbf{v}\mid\Theta_{i})$.

We parameterise the selection function $\mathbb{P}(\mathrm{S}\mid\mathbf{v},\Theta_{i})$ also as a bivariate GMM in colour-magnitude space with parameters $\widetilde{\boldsymbol{\boldsymbol{\epsilon}}_{i}}$

[TABLE]

The log likelihood is then given by

[TABLE]

As the selection function is a probability distribution, it must fall in the range $g_{i}^{\mathrm{SF}}(\mathbf{v}\mid\widetilde{\boldsymbol{\epsilon}_{i}},\Theta_{i})\in[0,1]$ for all $\mathbf{v}$. We reparametrise the selection function component weights as $\widetilde{\kappa_{ik}}=\mathrm{logit}\left(\frac{\widetilde{w_{ik}}}{\sqrt{|2\pi\widetilde{\boldsymbol{\Sigma}_{ik}}}}\right)$ with a uniform prior $\widetilde{\kappa_{ik}}\sim\mathrm{U}[-\inf,\inf]$. This constrains the component weights to the range $\widetilde{w_{ik}}\in\left[0,\sqrt{|2\pi\widetilde{\boldsymbol{\Sigma}_{ik}}|}\right]$ such that the maxima of any Gaussian component is less than or equal to one. The sum of selection function Gaussian mixture components must also be less than or equal to one everywhere. Since the maxima of a GMM is non-analytic, we find the roots of the gradient of the GMM using the hybr method from MINPACK-1 (More et al., 1980) implemented in scipy initialising at the mean of each GMM component. If the value of the GMM at any root is greater than one, the posterior probability is set to zero.

An NIW prior is used for $\widetilde{\boldsymbol{\mu}_{ik}},\widetilde{\boldsymbol{\Sigma}_{ik}}$ of the selection function.

We fit simultaneously for both the selection function parameters, $\widetilde{\boldsymbol{\epsilon}_{i}}$ and photometric density parameters, $\boldsymbol{\epsilon}_{i}$ where the prior distributions on $\boldsymbol{\epsilon}_{i}$ are the posteriors of the fit to the photometric sample in Section 2.2.1.

Combining Equation (19) with the priors on all parameters the posterior is given by

[TABLE]

To fit the selection function, we need to find $\mathrm{argmin}_{\boldsymbol{\epsilon}_{i},\widetilde{\boldsymbol{\epsilon}}_{i}}\left(-\ln\mathbb{P}\right)$. All terms in the posterior probability are analytically differentiable so we employ the ‘Truncated Newton’ method in scipy which uses the gradient of the probability distribution to converge on the minima. Parameters of the photometric density, $\boldsymbol{\epsilon}_{i}$ are initialised at the mean values of the prior distribution. The selection function parameters are initialised using two methods:

• •

$k$-means clustering with each spectroscopic sample star providing a weighted contribution of $\frac{1}{n_{i}^{\mathrm{phot}}(\mathbf{v}_{j}\mid\boldsymbol{\epsilon}_{i},\Theta_{i})}$. The reweighting is akin to approximating the selection function as $n^{\mathrm{spec}}/n^{\mathrm{phot}}$.

• •

Fit a GMM directly to the spectroscopic data using BayesianGaussianMixture with least informative Dirichlet and NIW priors.

The parameters are optimised for both initialisations and the one which leads to the largest posterior probability is used as the best fit. This helps to avoid some local optima in the posterior distribution.

In order to determine the optimal number of Gaussian components to use for the selection function, the algoritm is run for $\widetilde{K}=[1,K]$ components. We do not allow more components than the photometric density GMM as this would likely result in overfitting of the selection function. We use the number of components which minimizes the BIC

[TABLE]

where the number of degrees of freedom $d=(K+\widetilde{K})\times 6$ and $\mathcal{L}^{\mathrm{spec}}_{\mathrm{max}}$ is the likelihood corresponding to the best fit parameters. To generate a full posterior on the parameters of the photometric density and selection function, we run a set of $6\times(K+\widetilde{K})\times 2$ chains for 2000 iterations with emcee (Foreman-Mackey et al., 2013) initialising from a small ball around the best fit parameters. For every optima of the posterior distribution, there are a set of $\widetilde{K}!$ degenerate optima generated by reordering components. For finding the best fit parameters, it is unimportant which of the optima we sample around, however the posterior on the parameters should not be degenerate. We constrain the means of the photometric density and selection function components to maintain their respective orders in colour and apparent magnitude throughout the emcee iterations.

The models, their parameters and the prior distributions are summarised in Table 1 for reference. The choice of model hyperparameters is discussed further in Section 3.3.1.

2.3 Intrinsic coordinates

Chemodynamical models of the Milky Way make predictions for the metallicities, masses, and ages of stars at different positions. In order to test these models against observations, the selection function is required to account for the observation biases of the survey. To use the selection function, we require a transformation between observable coordinates (colour and apparent magnitude) and intrinsic coordinates (distance, metallicity, mass, and age). The selection function in terms of the intrinsic coordinates of the stars is

[TABLE]

where $s$ is the heliocentric distance, and $\mathrm{[M/H]}$, $\mathcal{M}_{\mathrm{ini}}$, and $\tau$ are the star’s metallicity, initial mass, and age respectively.

Here we follow a similar method for transforming between observable and intrinsic coordinate systems as described in Sanders & Binney (2015) and Das & Binney (2016). Any combination of $\mathrm{[M/H]}$, $\mathcal{M}_{\mathrm{ini}}$, and $\tau$ maps to a single set of coordinates, $(c,M)$, which are the colour and absolute magnitude of the star. Any values of $M$ and $s$ uniquely define $m$, the apparent magnitude of the star. Therefore any intrinsic coordinates, $(s,\mathrm{[M/H]},\mathcal{M}_{\mathrm{ini}},\tau)$ map to observable coordinates, $(c,m)$ and to a single value of the selection function, $\mathbb{P}(\mathrm{S}\mid\mathbf{v},\boldsymbol{\theta})$. The selection function in intrinsic coordinates for a star on field $i$ is therefore given by

[TABLE]

where $\mathbf{v}=(c,m)$. The maps are given by

[TABLE]

and

[TABLE]

where $F^{\mathrm{iso}}$ is the mapping introduced by the isochrones.

There are a variety of methods for generating the $F^{\mathrm{iso}}$ map. The most common method is to adopt the nearest isochrone to the values of age and metallicity provided. For this work, we linearly interpolate between isochrones, which improves the accuracy of the transformation.

We use the PARSEC isochrones (Bressan et al., 2012) on a grid of $353$ ages and $57$ metallicities. Ages in the range $-2.40\leq\log_{10}\left(\tau/\mathrm{Gyr}\right)\leq 1.12$ with a spacing of $0.01$ dex, and metallicities in the range $-2.192\leq\mathrm{[M/H]}\leq 0.696$ with a spacing of $0.051$ dex, were considered.

Every isochrone has a maximum initial mass above which a star of the given age and metallicity cannot exist. We generate a scaled initial mass coordinate so that each isochrone varies from $\widetilde{\mathcal{M}}_{\mathrm{ini}}\in[0,1]$,

[TABLE]

where $\mathcal{M}_{\mathrm{ini}}^{\mathrm{max}}(\tau,\mathrm{[M/H]})$ is the maximum initial mass value for a star of a given age and metallicity which we determined by linearly interpolating between maximum initial mass of isochrones as a function of $\tau,\mathrm{[M/H]}$. Likewise for $\mathcal{M}_{\mathrm{ini}}^{\mathrm{min}}(\tau,\mathrm{[M/H]})$.

The scaled mass is linearly interpolated along each isochrone. A single set of $\widetilde{\mathcal{M}}_{\mathrm{ini}}$ values is drawn for all isochrones which samples most heavily where the curvature of the isochrone in colour-apparent magnitude space is greatest. This parameterisation enables us to interpolate colour and absolute magnitude on a regular $\tau-[\mathrm{M/H}]-\widetilde{\mathcal{M}}_{\mathrm{ini}}$ grid. The interpolation is the mapping $F^{\mathrm{iso}}(\tau,\mathrm{[M/H]},\widetilde{\mathcal{M}}_{\mathrm{ini}}(\mathcal{M}_{\mathrm{ini}}))$.

3 Mock tests

To test the performance of the method presented in Section 2, we apply it to mock samples with known selection functions.

3.1 Galaxia mock catalogue

We generate our mock catalogue using the Galaxia code (Sharma et al., 2011). The thin disk is generated from the analytic Besançon model (Robin et al., 2003). Non-circular motions in the plane of the disk are introduced through the Shu (1969) DF, and the Bullock & Johnston (2005) N-body models simulate any substructure in the halo. The code synthesises a population of stars with coordinates, trajectories and intrinsic properties including age, metallicity and mass. Given the age, metallicity, initial mass and distance, Galaxia also calculates the colour and apparent magnitude of stars as observed from the Sun by employing the nearest Padova isochrone (Marigo et al., 2008) from a grid of 182 ages and 34 metallicities. The PARSEC (Bressan et al., 2012) isochrones however present a significant update on the Padova isochrones in terms of revisions to major input physics such as the equation of state, opacities, nuclear reaction rates, and inclusion of the pre-main sequence phase. The nearest isochrone calculation method is also less accurate than the interpolation method we employ with PARSEC isochrones discussed in Section 2.3. For these reasons, we recalculate the Galaxia apparent magnitudes and colours using our method.

We sample from Galaxia with an $H$-band magnitude limit, $m_{\mathrm{H}}<15$, similar to the limitations of the 2MASS survey (Skrutskie et al., 2006). The $H$-band magnitude limit was placed using Galaxia such that stars were filtered out based on the built-in Padova magnitudes. The difference between this and our own magnitude calculation leads to some stars dimmer than the magnitude cut being included and some brighter stars being excluded. This does not significantly affect our tests.

We will refer to this catalogue as the photometric sample.

3.2 Imposed selection function

The selection function is applied across two fields with coordinates

[TABLE]

Each field has a half-opening angle of $2\mathrm{deg}$ and hence a solid angle of $12.6\mathrm{deg}^{2}$. The locations of the fields in galactic coordinates are given in Figure 2 where the red and green shaded regions are fields 1 and 2 respectively and the greyscale background shows the number density of stars in the photometric sample.

For each field we apply three different selection functions:

• •

Flat: Selection boundaries of $m_{H}<13.5$ and $J-K>0.5$. The value of the selection function is $0.1$ within the boundaries and [math] outside.

• •

tanh:

[TABLE]

This is of a similar form to the selection function proposed for SEGUE G-dwarfs by Bovy et al. (2012).

• •

RAVE: We take 10 randomly selected fields from the RAVE survey (Kordopatis et al., 2013) along with any 2MASS (Skrutskie et al., 2006) stars on the same fields, bin the stars in $m_{H,\mathrm{2MASS}}-(J-K)_{\mathrm{2MASS}}$ and use the ratio of RAVE to 2MASS stars in each bin as the selection function probability.

The three applied selection functions are plotted in top panels of Figure 3. The red and green scatter points in Figure 2 are the stars selected by the Flat selection function to demonstrate the setup.

In the region of overlap of the fields, stars may be selected by either or both of the fields. We include every star selected by at least one field however we also record which fields each star was selected by as this is important for calculating the selection function. For instance, when deriving the selection function for field 1, we must use all stars which were selected by field 1 even if they were also selected by field 2. Likewise for field 2, we must include all stars selected by field 2 even if they were also selected by field 1. The effective double counting of selection functions in overlapping regions is accounted for by the union calculation as described in Section 2.1.1.

3.3 Results

The Galaxia photometric and spectroscopic catalogues are used to calculate the selection functions in observable coordinates, ${\mathbb{P}(\mathrm{S}\mid\boldsymbol{\theta},c,m)}$, and intrinsic coordinates, $\mathbb{P}(\mathrm{S}\mid\boldsymbol{\theta},s,\mathrm{[M/H]},\mathcal{M}_{\mathrm{ini}},\tau)$ following the methods described in Section 2.

We first discuss the values of hyperparameters used for the GMM priors. We will present and test the posterior for field 1 as a function of colour and apparent magnitude. We then show the fit to both fields in observable coordinates and position on the sky which tests how the model handles overlapping fields. Finally we present the selection function in intrinsic coordinates.

3.3.1 Prior hyperparameters

For the prior distributions on the photometric density and selection function we require hyperparameter values which encode our prior knowledge about the distribution. We assume that we know all data falls within the apparent magnitude and colour ranges $m_{H}\in[4,15]$ and $J-K\in[-0.1,1.2]$.

For the NIW hyperparameters on the photometric density fit we use $\mathbf{m}_{0}=\mathrm{med}\left[\mathbf{v}_{\mathrm{min}},\mathbf{v}_{\mathrm{max}}\right]=(9.5,0.55)$ which is at the center of our prior range. We set $\lambda=10^{-4}$. By choosing a small value for $\lambda$ we ensure that the prior is extremely uninformative on the mean. For the covariance we choose $\boldsymbol{\Psi}=\mathrm{diag}((\frac{\mathbf{v}_{\mathrm{max}}-\mathbf{v}_{\mathrm{min}}}{2})/5)^{2})=\mathrm{diag}(1.21,0.0169)$ such that the standard deviation is $1/5$ the prior range in $m_{H}$, $J-K$. The precise choice of $1/5$ is somewhat arbitrary and chosen as it leads to reasonable numbers of Gaussian components in the photometric density and selection function GMMs in all cases. The degrees of freedom is chosen as $\nu=2$ which is the most uninformative choice. For the Dirichlet prior we choose a concentration prior of $\alpha_{k}=1/K$ for all components which is the least informative prior.

We use the same hyperparameters for the NIW prior on the selection function as those used for the photometric density fit. This is because we have no more information about the selection function than we did about the photometric density other than that it is only important in the given colour and magnitude range. Having found the best fit parameters, we also run a set of emcee chains as discussed in Section 2.2.2 to generate posterior samples on all of the photometric density and selection function parameters.

The choices of hyperparameter values are summarised in the right-hand column of Table 1.

3.3.2 Observable coordinates

In this section we demonstrate the selection function fit to field 1 of the mock sample in colour-apparent magnitude space. To minimize the BIC, the optimal number of photometric density GMM components was 14 whilst the Flat, tanh and RAVE selection functions were fit with 2, 2 and 1 components.

For the RAVE mock, the distribution of photometric and spectroscopic points and the GMM distribution fits are given in Figure 4. For a successful fit we would expect the centers of $68\%$ of histogram bins to fall within the Poisson uncertainties of the model and we see that this appears to be the case for our fit to the RAVE selection function against both colour and apparent magnitude. The same number density distributions for all three mock selection functions are shown in Figure 5. For the Flat model (top panel) the smooth GMM selection function fails to correctly fit the sharp cuts at $m_{H}=13.5$ and $J-K=0.5$. This is a predictable limitation of attempting to fit a discontinuous model with a smooth function. In all other cases the histograms of the spectroscopic sample mostly fall within the Poisson noise uncertainties of the products between our fits to the photometric density and selection function. This provides a qualitative demonstration that the method is providing reasonable results.

For a more quantitative assessment we use the Kolmogorov-Smirnov (KS) test. We randomly draw 1000 samples from the last 500 iterations of each emcee fit to the spectroscopic sample. The KS probability (p-value of the KS test) is the probability that the spectroscopic sample is consistent with being drawn from the given probability distribution, in this case the product of the photometric and selection function GMMs which is itself a GMM with $K\times\widetilde{K}$ components.

The KS statistics are computed in 1D as a function of colour and apparent magnitude separately. For a good fit to the data, the KS probabilities will be uniformly distributed in the range $[0,1]$. If the probability is skewed towards [math], the method has underfit or failed to fit the data. If the probability is skewed towards $1$, the method has overfit the data and hence the selection function would not generalise well to a new dataset. The most extreme example of this would be placing delta functions on every star which would achieve a KS probability of 1.

We expect that the KS statistic should demonstrate some overfitting since the test does not consider the prior information provided by the photometric sample.

Histograms and cumulative distributions of the KS tests for field 1 against colour and apparent magnitude are shown in Figure 6. Due to the sharp cut off in the Flat model, the KS probabilities demonstrate that the GMM underfits the selection function as a function of magnitude. Curiously, the same problem is not seen as a function of colour where the Flat model is almost perfectly fit. KS tests are only weakly sensitive to the wings of the distribution and we postulate that this is why the sharp cutoff issue is not picked up in colour space. The RAVE and tanh models provide reasonable fits as a function of apparent magnitude however the RAVE model is heavily overfit as a function of colour. The cause of this overfitting is unclear.

To test the normalisation of the photometric density and selection function we compare the integral over colour and apparent magnitude of the GMMs with the number of objects in the spectroscopic sample. For large samples, we would expect the difference between the integral and true count to be Gaussian distributed with variance equal to the size of the sample. Figure 7 shows the distribution of these integrals from the 1000 emcee parameter draws re-centered and re-scaled by the spectroscopic sample number counts. All distributions are centered around zero which implies that the fits are well normalised. The fits to the RAVE model show a tight distribution which is another indication of overfitting to the spectroscopic data.

3.3.3 Overlapping fields

On field 2 the photometric density is fitted with 14 GMM components and the Flat, tanh and RAVE selection functions with 2, 3 and 1 components respectively to minimize the BIC. We will now only consider the selection function resulting from the combination of fields 1 and 2.

The selection function is now a superposition of two GMMs as a function of galactic coordinates. As such the KS one-sample statistic can no longer be used as we don’t have an analytic distribution to fit. Instead we draw random samples from the photometric catalogue with the probability of inclusion given by the model selection function as a function of galactic coordinates, colour and apparent magnitude. We compute the two-sample KS probability against the spectroscopic catalogue which gives the probability that samples are drawn from the same probability distribution.

In Figure 8 we show the distribution of KS probabilities for fields 1 and 2 with all combinations of the Flat, tanh and RAVE selection functions applied. The Flat+Flat model (blue solid line) systematically underfits the data in both apparent magnitude and colour which, as discussed earlier, is caused by the sharp cuts in the selection function which a smooth GMM cannot correctly reproduce. All other combinations of fields produce close to uniform distributions of KS probabilities with weak under and overfitting against different coordinates.

To test the normalisation of the selection function we compare the sum of the selection function over the stars in the photometric sample with the number in the spectroscopic sample. This sum is the mean sample size which would be produced by drawing many selection function weighted samples from the photometric catalogue. The re-centered and re-scaled distributions for all field combinations are shown in Figure 9. We see significantly less overfitting here with all distributions reproducing the expected Poisson noise. Only the Flat+tanh model (blue dashed) systematically overestimates the normalisation however the offset is only at the one standard deviation level so we don’t consider this to be a significant issue.

These tests demonstrate the effectiveness of the union method described in Section 2.1.1 for evaluating the selection function of overlapping fields of the spectroscopic survey. This also shows that the method provides good fits to selection functions with very different properties. A caveat is that the GMM is not perfectly suited to fitting selection functions with discontinuous changes however even in these cases reasonable fits can still be achieved.

We also show here the power of our method for calculating selection functions of combined surveys. For example we could apply this method to generate a single selection function for the combined APOGEE and RAVE catalogues given their individual selection functions.

3.3.4 Intrinsic coordinates

Our final test is on selection function as a function of intrinsic coordinates (age, metallicity, mass and distance) using the mapping laid out in Section 2.3. We test the selection function for the two field sample with the tanh and RAVE selection functions applied to fields 1 and 2 respectively. A histogram for the model fit is generated by taking the binned sum of the selection function probabilities of the photometric sample and is shown by the orange dashed histogram in Figure 10. The blue solid line shows the histogram for the spectroscopic sample and this falls within Poisson noise uncertainties of the model in the majority of bins.

Similar to Section 3.3.3, we generate 1000 mock samples from the photometric catalogue weighted by the inferred selection function and calculate the goodness of fit from the two-sample KS statistic. The distribution of KS probabilities is given in Figure 11. As we saw in observable coordiantes, the distribution of KS probabilities is near uniform with very weak under fitting against distance and overfitting in initial mass and metallicity. This demonstrates that our method is also well suited to determining selection functions in intrinsic coordinates.

4 Discussion

In this paper, we have introduced a novel method for deriving selection functions of spectroscopic surveys. We have also demonstrated the success of this method on a set of test cases using Galaxia. In this section, we discuss some key points to consider when applying the method, and detail potential improvements to come.

4.1 Choice of photometric catalogue

When calculating the selection function of any spectroscopic survey, a photometric survey needs to be specified which may be assumed complete in the region of the colour-apparent magnitude space explored. The choice of photometric catalogue is dependent on the characteristics of the survey. The photometric catalogue should cover the whole footprint of the spectrograph or else a combination of photometric catalogues should be used. It is beneficial to use a photometric catalogue with observing bands closely matching the spectrograph wavelength range as this enables prior information on the spectrograph’s selection limitations to be applied more easily.

That said, particularly for low-latitude observations, dust attenuation is a significant factor, which suggests that infrared photometric surveys may be more appropriate. We will discuss the inclusion of dust attenuation to the intrinsic selection function in Section 4.3.

Gaia DR2 (Gaia Collaboration et al., 2018) represents the largest survey of the Milky Way to date and will be a complete photometric survey for the magnitude ranges of many spectrographs, particularly in high latitude fields. The selection of Gaia DR2 as a function of $l,b,m_{G}$ will be provided by Boubert & Everall (2020a, in prep) and can be used to test whether it is complete in the required magnitude range.

4.2 Error convolution

The selection function we have derived here is the probability of selection given measured properties of the stars. By convolving the likelihood with the measurement uncertainty of the photometries of each star, we can in principle derive the selection function given true properties of the stars. By virtue of defining the selection function as a GMM in colour-apparent magnitude space the convolution is analytic and as such the calculation is computationally feasible.

We do not present this here, but consider it as a potential avenue to pursue in the future.

4.3 Dust attenuation

In this work we have not considered the impact of interstellar dust on the derived selection function. The observable selection function is dependent on dust attenuation only through it’s effect on colour and apparent magnitude. Our inferred selection function given colour and apparent magnitude will be unchanged.

Dust attenuation changes the mapping from intrinsic coordinates to observable coordinates discussed in Section 2.3. The effect of dust on the transformation from absolute to apparent magnitude of the star is given by

[TABLE]

where $x$ represents the observation band being used for apparent magnitude in the selection function. Likewise the colour will also need to be corrected for dust reddening in the mapping from intrinsic coordinates to observables.

This will be an especially important consideration for low-latitude fields. To do this, we require an adequate 3D all-sky dust map.

Bovy et al. (2016) construct a composite map that patches together maps of Marshall et al. (2006), Sale & Magorrian (2014) and Green et al. (2015). A significant improvement on this is provided by Green et al. (2019) for $\mathrm{dec}>-30\mathrm{deg}$ which is extrapolated to the whole sky using Gaia photometries in Boubert & Everall (2020b, in prep).

We can include these maps in our model as an extinction and reddening term in the intrinsic to observable coordinate mapping described in Section 2.3.

4.4 Gaia RVS

We have not considered the Gaia Radial Velocity Spectrograph (RVS Katz et al., 2019) in this paper as the strong dependence on the Gaia scanning law means that our method is not the best for this catalogue. In Boubert & Everall (2020b, in prep) we develop an independent method for evaluating the selection function of Gaia RVS which properly treats the intricate behaviour of the scanning law. These works will prove complimentary to one another for studies of the chemo-dynamic properties of the Milky Way.

4.5 Future work

We will be applying this method to a series of currently available multifibre spectroscopic datasets including APOGEE DR14 (Abolfathi et al., 2018), LAMOST DR5 (Zhao et al., 2012), RAVE DR5 (Kunder et al., 2017), Gaia-ESO DR3 (Gilmore et al., 2012), GALAH DR2 (Buder et al. (2018)) and SEGUE (Yanny et al., 2009) and their future data releases. Further into the future, we will apply this algorithm to the first data releases of upcoming multifibre spectrographs, WEAVE (Dalton et al., 2012), MOONS (Cirasuolo et al., 2014) and 4MOST (de Jong et al., 2012).

All new selection functions will be made publicly available on the Python module selectionfunctions (https://github.com/DouglasBoubert/selectionfunctions) produced by Boubert & Everall (2020a, in prep) based on the dustmaps source code (Green, 2018).

5 Conclusions

We have developed a Bayesian model for empirically determining the selection function of multi-fibre spectrographs, where there exists a complete photometric survey in the same region of observable (colour-apparent magnitude) parameter space.

The method improves on previous works by modelling the selection function with a Gaussian mixture model that is fit to the data through a Poisson likelihood function. This generates a selection function which accounts for Poisson noise in low-count data. This approach also allows us to define the uncertainties in our selection function by analysing the posterior distribution on the model parameters.

We further incorporate a union calculation which allows the selection function to be calculated in regions of sky where fields partially or fully overlap. This can be applied to merged catalogues of independent surveys to produce a combined selection function. In an era where large amounts of spectroscopic data are becoming available from many independent observatories, each with their own observational limitations, combining surveys can hugely enhance our understanding of the Milky Way. We can also apply our method to subsamples of any catalogue if analysis is being done on a more specific or constrained stellar population.

Finally, we present a method of translating selection functions from observable colour-apparent magnitude coordinates into intrinsic coordinates of age, metallicity, mass and distance using the PARSEC isochrones (Bressan et al., 2012). This allows a deeper insight into the effects of the selection function on stellar parameters.

We have demonstrated the effectiveness of our method on a mock catalogue generated using the Galaxia (Sharma et al., 2011) population synthesis code. We are successful in reproducing the applied selection function within Poisson noise uncertainty. Using KS tests in one dimension we show that our method produces good fits to the data only struggling where the model selection functions undergo large discontinuous transitions.

Our code is made publicly available as a Python repository at https://github.com/AndrewEverall/seestar.git which will continue to be improved and updated with new developments. We will apply our method to several spectroscopic catalogues and make these selection functions publicly available.

Acknowledgements

AE is grateful to the Galactic Dynamics group in Oxford for providing the support and funding necessary to develop this work. AE in particular thanks Prof James Binney, Dr John Magorrian, and Dr Ralph Schönrich for regular discussions, which led to many improvements. AE thanks Zephyr Penoyre and Stephen Thorp for useful discussions on optimisation techniques. We are extremely greatful to the referee who provided detailed comments which lead to significant improvements to the paper and the method.

Appendix A Union expansion

Considering the simplest case of two overlapping fields, A and B, Equation (5) expands to

[TABLE]

In most surveys we are considering, there will be substantially more than two fields. For $M$ fields, we expand the union using inclusion-exclusion principle.

[TABLE]

where we have dropped the conditionals such that $\mathbb{P}(\mathrm{S}_{i})\equiv\mathbb{P}(\mathrm{S}_{i}\mid\boldsymbol{\theta},\mathbf{v})$ for ease of notation.

Assuming independence of observations from different fields (i.e. the event that a star is selected on field A is independent of whether it has been observed on field B) we can expand the joint probability as the product of the probability of each event.

[TABLE]

Using Equation 4 we expand the conditional selection probabilities in terms of the event that the positional coordinates lie on the field

[TABLE]

where $\mathbb{P}(\Theta_{i}\mid\boldsymbol{\theta})=1$ if $\boldsymbol{\theta}$ is on field i and $\mathbb{P}(\Theta_{i}\mid\boldsymbol{\theta})=0$ otherwise. Therefore any joint probability terms with fields which don’t contain $\boldsymbol{\theta}$ will vanish from Equation LABEL:eq:expansion.

We can simplify Equation (LABEL:eq:expansion) for some specific circumstances:

1. (1)

$\boldsymbol{\theta}$ is not located on any fields. $\mathbb{P}(\Theta_{i}\mid\boldsymbol{\theta})=0\quad\forall\quad i$. All terms in the expansion are [math]

[TABLE] 2. (2)

$\boldsymbol{\theta}$ is located on only one patch, $A$. $\mathbb{P}(\Theta_{i}\mid\boldsymbol{\theta})=0\quad\forall\quad i\neq A$, $\mathbb{P}(\Theta_{A}\mid\boldsymbol{\theta})=1$

[TABLE] 3. (3)

$\boldsymbol{\theta}$ is on the intersection between two fields denoted by $A$ and $B$:

[TABLE]

Appendix B Poisson likelihood

We start from the likelihood of observing a particular object with coordinates $v_{i}$ given that we only pick one point from the density function, $\lambda$

[TABLE]

Expanding this to observations of a population of N objects

[TABLE]

The likelihood of the data is then given by

[TABLE]

The denominator here is independent of the parameters of the density model so the likelihood of our dataset, $\{v\}$ is given by

[TABLE]

And taking the log likelihood

[TABLE]

Appendix C Model Structure.

In Table 2 we provide short descriptions of the notation followed in this paper to help with following the method.

In Figure 12 we present a Bayesian network (or acyclic directed graph) of the method for calculating the selection functions for each field as a function of colour and apparent magnitude. Red boxes contain parameters and hyperparameters of the model. Green ellipses are conditional probabilities from which parameters and data are sampled. The blue double circles are the observable data. In summary, the method consists of maximising the product off all probability distributions in green ellipses. The best fit hyperparameters then define the posterior parameters of the photometric density and selection function GMMs. In our method this is achieved as a two stage process, first fitting the photometric data to determine posterior GMM parameters for the photometric density. Subsequently the photometric GMM parameter posteriors are provided as priors to the spectroscopic likelihood function.