seestar: Selection functions for spectroscopic surveys of the Milky Way
Andrew Everall Payel Das


firstpage–lastpage
Multi-variant scheduling of critical time-triggered communication in incremental development process: Application to FlexRay
    Jan Dvořák12,
           Zdeněk Hanzálek2

        1Department of Control Engineering, Faculty of Electrical Engineering, CTU in Prague

        2Industrial Informatics Research Center, Czech Institute of Informatics, Robotics and Cybernetics, CTU in Prague
    Received: 28 September 2018 / Accepted: 5 February 2019
=====================================================================================================================================================================================================================================================================





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.





surveys – instrumentation: spectrographs –Galaxy: stellar content – methods: statistical










§ 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 <cit.> 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. <cit.> and <cit.> 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 <cit.>. Similar methods are used in determining the selection function for the Radial Velocity Experiment <cit.> by <cit.> and for the Large Sky Area Multi-Object Fibre Spectroscopic Telescope <cit.> Spectroscopic Survey of the Galactic Anticentre <cit.> by <cit.>. <cit.> perform a similar analysis of the Gaia-ESO survey <cit.>, taking into account the survey's observing strategy.

<cit.> generate a combined selection function for LAMOST <cit.>, RAVE <cit.> and the Tycho-Gaia Astrometric Solution (TGAS, ) by binning the fields in colour, apparent magnitude, and distance using a synthetic galaxy catalogue as the assumed complete sample. <cit.> and <cit.> consider the dependence of the selection function on apparent magnitude for G-type dwarfs in the Sloan Extension for Galactic Understanding and Exploration <cit.> survey and the Apache Point Observatory Galaxy Evolution Experiment <cit.>. 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. <cit.> also present a selection function in distance that takes into account the dust extinction along the line of sight.

Finally, <cit.> and <cit.> 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, ).

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 <cit.>, <cit.>, and <cit.>. With the exception of <cit.>, no methods are presented for combining selection functions.


In this paper, we build on previous work to create [https://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 <cit.>.

In Section <ref>, 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 <ref>. Finally, we discuss our results and potential future developments and applications of this method in Section <ref>. Section <ref> 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 <ref> provides a table of parameter descriptions for reference along with a graphical Bayesian network to help with following our notation and understanding the method.


 

§ METHOD


The selection function (|𝐱) of a stellar survey is the probability of a star being in the survey, , given the star's coordinates 𝐱. 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


    (|𝐱) = (𝐱|)   () /( 𝐱 ).




(𝐱) d^n x ≡ f(𝐱) d^n x is the probability that a star chosen at random has coordinates within the volume element, d^nx, where f(𝐱) is the distribution function of stars in the Milky Way. (𝐱|)  d^nx ≡ f^spec(𝐱|)  d^nx is the probability that a star picked at random from the survey has coordinates in d^nx, where f^spec(𝐱) is the distribution function of observed stars. Finally, () = 𝔼[N^spec] / 𝔼[N] 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^phot(𝐱) = f(𝐱) and N^phot=N and substitute into Equation (<ref>).


    (|𝐱) =  f^spec(𝐱) 𝔼[N^spec]/ f^phot(𝐱) 𝔼[N^phot] =
    𝔼[n^spec(𝐱)]/𝔼[n^phot(𝐱)],


where 𝔼[n^spec(𝐱)] and 𝔼[n^phot (𝐱)] are the expected number densities of stars in the spectroscopic and photometric surveys respectively, at coordinates 𝐱.

The observation of a star is dependent on the star's observable properties, 𝐱 = (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, = (l, b) and photometric properties, 𝐯=(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 <ref>. On each patch, the selection function is calculated as a function of 𝐯 as described in Section <ref>. In Section <ref>, the observed coordinates are transformed into intrinsic coordinates, 𝐯(s, , ℳ_ini, τ), where s, , ℳ_ini, τ are distance, metallicity, initial mass, and age respectively. These coordinates underpin the description of the chemodynamical distribution of stars in the Milky Way.




 §.§ 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 Θ_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


    (Θ_i  |) =
    
    1    if  is located on patch i,
    
    0    otherwise.






  §.§.§ Defining a patch



We assume that the selection function given 𝐯 is independent of  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, S/N∼√(λ) where λ 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 <ref> 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, 𝐯. If a coordinate does not lie in a field, the selection function is 0.

The selection function for field i is given by


    (_i|, 𝐯) = (_i |Θ_i, 𝐯)(Θ_i |)


where _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 <ref>, 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.


    (|, 𝐯 )  = (⋃_i=1^M_i | , 𝐯)


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


    (⋃_i^M _i | , 𝐯) = ∑_k=1^M (-1)^k+1[   ∑_1≤ i_1 < ... < i_k ≤ M(_i_1, _i_2 ...i_k) ]

 
where 


    (_1, _2, ...) = ∏_i = 1, 2, ...(_i |Θ_i, 𝐯) (Θ_i |).


assuming the events _1, _2... are independent. Appendix <ref> 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 (<ref>) 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 𝐯. These are calculated for a given field, i. Having calculated the selection function for each field as a function of 𝐯, they are combined using the method above to achieve the full selection function for the entire survey.




 §.§ Selection function in observable coordinates



The probability of a star being selected by field i given 𝐯, and given that the star lies on the field (i.e. (Θ_i|) = 1) is


    (_i|𝐯,Θ_i) = (𝐯|_i, Θ_i)   (_i|Θ_i)/(𝐯|Θ_i) .


(𝐯|Θ_i) d𝐯≡ f (𝐯|Θ_i) d𝐯 is the probability that a star chosen at random on field i has coordinates within the d𝐯 volume element, where f(𝐯|Θ_i) is the distribution of stars in the Milky Way inside the cone projecting onto field i.

(𝐯|_i, Θ_i) d𝐯≡ f^spec(𝐯|Θ_i)  d𝐯 is the probability that a star observed on field i of the spectroscopic survey has coordinates within the d𝐯 volume element. Finally, (_i) = 𝔼[N_i^spec] / 𝔼[N_i] 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^phot(𝐯 | Θ_i) = f(𝐯 | Θ_i) and N_i^phot = N_i. Substituting into Equation (<ref>)


    (_i|Θ_i,𝐯) = f_i^spec (𝐯|Θ_i)𝔼[N_i^spec]/f_i^phot(𝐯|Θ_i) 𝔼[N_i^phot]
    = 𝔼[n_i^spec(𝐯|Θ_i)]/𝔼[n_i^phot(𝐯|Θ_i)],


where n_i^spec(𝐯|Θ_i)=f_i^spec (𝐯|Θ_i)N_i^spec is the number density of stars observed by the survey on field i and n_i^phot(𝐯|Θ_i)=f_i^phot (𝐯|Θ_i)N_i^phot is the number density of stars in the Milky Way on the cone projecting onto field i.





  §.§.§ Number density of photometric sample



We start by calculating the expected number density of stars in the Milky Way on field i, 𝔼[n_i^phot(𝐯|Θ_i)]. The choice of photometric survey is discussed in Section <ref>.


The stars in the photometric survey represent a Poisson realisation of the smooth underlying number density function, n^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


    G(𝐯|μ, Σ) = 1/√(|2πΣ|)exp(- (𝐯-μ)^T  Σ^-1  (𝐯-μ)/2),


where μ∈ℝ^2 is the mean of the bivariate Gaussian in colour and magnitude and Σ∈ℝ^2×2 is the symmetric covariance matrix. Using this, we can write


    n_i^phot(𝐯|_i,Θ_i) = ∑_k=1^K w_i k   G(𝐯|μ_i k, Σ_i k),


where _i ∈ℝ^K× 6 for the K components of the GMM, each with six parameters defining the bivariate Gaussian component w, μ^0, μ^1, Σ^00, Σ^11 and Σ^01 (Σ is symmetric so Σ^10=Σ^01).

We then find the parameters _i by maximising the likelihood of the photometric catalogue stars. The Poisson likelihood is derived from Poisson count probabilities to give

    ln(ℒ^phot(X_i^phot | _i)) = -∫   d𝐯  n_i^phot(𝐯|_i, Θ_i) 
       + ∑_j=1^N_i^photlog( n_i^phot(𝐯_j |_i, Θ_i) ),


where 𝐯_j are the colour and apparent magnitude of star j in the photometric catalogue on field i (see Appendix <ref>).

We reparameterise the weights as normalised `loadings', π_i k = w_i k/N where N = ∑_k=1^K [w_i k]. 

    ln (ℒ^phot(X_i^phot | _i)) ∝ -    N + N^phot_iln(N) 
       +∑_j=1^N^phot_iln(n̂^phot_i (𝐯_j |_i, Θ_i)),

where n̂^phot_i = n^phot_i/N. The normalisation, N is now a parameter of the model whilst the loadings provide K-1 free parameters under the constraint ∑_k=1^K π_i k = 1.

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


    NIW(   μ, Σ | 𝐦_0, λ, Ψ, ν) 
    ∝     𝒩(μ|𝐦_0, 1/λΣ)  
    |Σ|^-(ν + p + 1)/2exp[-1/2(ΨΣ^-1)],

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 μ and Σ will also be represented by a NIW distribution. The chosen hyperparameters (𝐦_0, λ, Ψ, ν) are dependent on the prior information about the data. Choosing small values of λ and setting ν = p provides a least informative prior on the parameters. We discuss a good choice of hyperparameters in our mock tests in Section <ref>.
 
The loadings are evaluated with a Dirichlet prior,

    DIR(π_1, π_2...π_K | α_1, α_2...α_K) ∝∏_k=1^Kπ_k^α_k - 1

which satisfies the constraint ∑_k=1^K [π_i k] = 1.
The Dirichlet distribution is the conjugate prior for membership probability of a categorical variable. Setting concentration hyperparameters α_i = α = 1/K  ∀  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 <cit.> 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 α_k', 𝐦_0k', λ_k', Ψ_k' and  ν_k' 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 (<ref>) where the uniform prior doesn't provide a contribution. This is proportional to a Poisson distribution with mean N^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 <cit.> for the best fit parameters (taken as the expected values of the posteriors)

    BIC =ln(N_i^phot)d - 2ln(ℒ^phot_max)

where d=K× 6 (the number of parameters in our model) and ℒ^phot_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.





  §.§.§ Selection function in observed coordinates





The selection function in observed coordinates for field i is (_i |𝐯, Θ_i). The stars in the spectroscopic catalogue represent a Poisson realisation of the product of the distribution function and the selection function


    n_i^spec(𝐯|Θ_i) = n_i^phot(𝐯|Θ_i)  × (_i |𝐯, Θ_i).


Therefore, we can once more use the Poisson likelihood in Equation (<ref>), replacing n^phot(𝐯|Θ_i) with n^spec(𝐯|Θ_i).

We parameterise the selection function (|𝐯, Θ_i) also as a bivariate GMM in colour-magnitude space with parameters _i 

    g_i^SF(𝐯|_i,Θ_i) = ∑_k=1^Kw_i k   G(𝐯|μ_i k, Σ_i k),

The log likelihood is then given by

    ln(ℒ^spec(X_i^spec |    _i, _i)) ∝
    
       -∫d𝐯     n_i^phot(𝐯|_i, Θ_i)  g_i^SF(𝐯|_i,Θ_i) 
    
       +   ∑_j=1^N_i^speclog( n_i^phot(𝐯_j |_i, Θ_i)  g_i^SF(𝐯_j |_i,Θ_i)   ).

  
As the selection function is a probability distribution, it must fall in the range g_i^SF(𝐯|_i,Θ_i) ∈ [0,1] for all 𝐯. We reparametrise the selection function component weights as κ_i k = logit(w_i k/√(|2πΣ_i k)) with a uniform prior κ_i k∼U[-inf, inf]. This constrains the component weights to the range w_i k∈[0, √(| 2πΣ_i k|)] 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 <cit.> 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 μ_i k, Σ_i k of the selection function.

We fit simultaneously for both the selection function parameters, _i and photometric density parameters, _i where the prior distributions on _i are the posteriors of the fit to the photometric sample in Section <ref>.

Combining Equation (<ref>) with the priors on all parameters the posterior is given by

    ln((_i, _i | X_i^spec)) = ln(   ℒ^spec(X_i^spec | _i, _i))
    
    -    N + N_i^photln(N) 
    
    +    DIR(π_1, π_2...π_K | α_1', α_2'...α_K') 
    
    +    ∑_k=1^K [NIW(μ_i k, Σ_i k | 𝐦_0ik', λ_ik', Ψ_ik', ν_ik') ]
     
    +    ∑_k=1^K[NIW(μ_ik, Σ_ik | 𝐦_0, λ, Ψ, ν)].

To fit the selection function, we need to find argmin__i, _i(-ln). 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, _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  1/n_i^phot(𝐯_j|_i,Θ_i). The reweighting is akin to approximating the selection function as n^spec/n^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 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

    BIC =ln(N_i^spec)d - 2ln(ℒ^spec_max)


where the number of degrees of freedom d=(K+K)×6 and ℒ^spec_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 × (K + K) × 2 chains for 2000 iterations with  <cit.> initialising from a small ball around the best fit parameters.
For every optima of the posterior distribution, there are a set of 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  iterations.

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




 §.§ 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


    (|  s, , _ini, τ, θ),


where s is the heliocentric distance, and , _ini, and τ 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 <cit.> and <cit.>. Any combination of , _ini, and τ 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, , _ini,τ) map to observable coordinates, (c, m) and to a single value of the selection function, (|𝐯, θ). The selection function in intrinsic coordinates for a star on field i is therefore given by


    (|  s, τ, , _ini, θ ) = (|𝐯(s, τ, , _ini), θ).


where 𝐯=(c,m). The maps are given by


    m = M + 5log_10(s/10 pc)


and

    (c, M) = F^iso(τ, , _ini),


where F^iso is the mapping introduced by the isochrones.

There are a variety of methods for generating the F^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 <cit.> on a grid of 353 ages and 57 metallicities. Ages in the range -2.40≤log_10(τ/Gyr) ≤ 1.12 with a spacing of 0.01 dex, and metallicities in the range -2.192 ≤≤ 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 _ini∈  [0, 1],


    _ini = _ini - _ini^min/_ini^max - _ini^min,


where _ini^max(τ,) 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 τ,. Likewise for _ini^min(τ,). 

The scaled mass is linearly interpolated along each isochrone. A single set of _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 τ - [M/H] - _ini grid. The interpolation is the mapping F^iso(τ, , _ini(_ini)).









§ MOCK TESTS






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






 §.§  mock catalogue





We generate our mock catalogue using the  code <cit.>. The thin disk is generated from the analytic Besançon model <cit.>. Non-circular motions in the plane of the disk are introduced through the <cit.> DF, and the <cit.> 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,  also calculates the colour and apparent magnitude of stars as observed from the Sun by employing the nearest Padova isochrone <cit.> from a grid of 182 ages and 34 metallicities. The PARSEC <cit.> 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 <ref>. For these reasons, we recalculate the  apparent magnitudes and colours using our method.

We sample from  with an H-band magnitude limit, m_H<15, similar to the limitations of the 2MASS survey <cit.>. The H-band magnitude limit was placed using  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.













 §.§ Imposed selection function


The selection function is applied across two fields with coordinates

    l, b =
    
    30, 60    (Field 1),
    
    30, 61    (Field 2).


Each field has a half-opening angle of 2deg and hence a solid angle of 12.6deg^2. The locations of the fields in galactic coordinates are given in Figure <ref> 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 0 outside.

  * tanh:
    
    (S|m_H, J-K) = 0.1    ×(1- tanh(m_H - 13.5/e^-2))/2 
       ×(1- tanh(J-K - 0.5/e^-2))/2.

This is of a similar form to the selection function proposed for SEGUE G-dwarfs by <cit.>.

  * RAVE: We take 10 randomly selected fields from the RAVE survey <cit.> along with any 2MASS  <cit.> stars on the same fields, bin the stars in m_H, 2MASS - (J-K)_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 <ref>. The red and green scatter points in Figure <ref> 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 <ref>.





 §.§ Results



The  photometric and spectroscopic catalogues are used to calculate the selection functions in observable coordinates, (|, c, m), and intrinsic coordinates,  following the methods described in Section <ref>. 

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.




  §.§.§ 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 ∈ [4, 15] and J-K ∈ [-0.1, 1.2].

For the NIW hyperparameters on the photometric density fit we use 𝐦_0 = med[𝐯_min, 𝐯_max] = (9.5, 0.55) which is at the center of our prior range. We set λ = 10^-4. By choosing a small value for λ we ensure that the prior is extremely uninformative on the mean. For the covariance we choose Ψ = diag((𝐯_max-𝐯_min/2)/5)^2) = 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 ν = 2 which is the most uninformative choice. For the Dirichlet prior we choose a concentration prior of α_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  chains as discussed in Section <ref> 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 <ref>.

 


  §.§.§ 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 <ref>. 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 <ref>. 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  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×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 0, 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 <ref>. 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 <ref> shows the distribution of these integrals from the 1000  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.


 


  §.§.§ 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 <ref> 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 <ref>. 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 <ref> 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.


 


  §.§.§ 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 <ref>. 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 <ref>. 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 <ref>, 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 <ref>. 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.






§ 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 . In this section, we discuss some key points to consider when applying the method, and detail potential improvements to come.



 §.§ 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 <ref>.

Gaia DR2 <cit.> 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.



 §.§ 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.



 §.§ 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 <ref>. The effect of dust on the transformation from absolute to apparent magnitude of the star is given by

    m_x = M_x + 5log_10(s/10 pc) + A_x(l, b, s),

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. 

<cit.> construct a composite map that patches together maps of <cit.>, <cit.> and <cit.>. A significant improvement on this is provided by <cit.> for dec>-30deg 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 <ref>.




 §.§ Gaia RVS


We have not considered the Gaia Radial Velocity Spectrograph <cit.> 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.



 §.§ Future work


We will be applying this method to a series of currently available multifibre spectroscopic datasets including APOGEE DR14 <cit.>, LAMOST DR5 <cit.>, RAVE DR5 <cit.>, Gaia-ESO DR3 <cit.>, GALAH DR2 (<cit.>) and SEGUE () and their future data releases. Further into the future, we will apply this algorithm to the first data releases of upcoming multifibre spectrographs, WEAVE <cit.>, MOONS <cit.> and 4MOST <cit.>.


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 <cit.>.





§ 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 <cit.>. 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  <cit.> 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.



mnras










§ UNION EXPANSION

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


    (|, 𝐯) = (_A|, 𝐯) +(_B|, 𝐯) - (_A, _B|, 𝐯).


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.


    (⋃_i=1^M _i ) =    ∑_i=1^M (_i) -∑_i=2^M∑_j=1^i-1(_i, _j)  
       +∑_i=3^M∑_j=2^i-1∑_k=1^j-1(_i,_j,_k) - ... 
        ...  + (-1)^M+1(_1, _2...._M) 
    
     =    ∑_k=1^M (-1)^k+1[   ∑_1≤ i_1 < ... < i_k ≤ M(_i_1, _i_2 ...i_k) ],

where we have dropped the conditionals such that (_i) ≡(_i |, 𝐯 ) 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.

    (_1, _2, ...|, 𝐯) = ∏_i = 1, 2, ...(_i |, 𝐯)


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

    (_1, _2, ...|, 𝐯) = ∏_i = 1, 2, ...(_i |Θ_i, 𝐯) (Θ_i |)

where (Θ_i |) =1 if  is on field i and (Θ_i   |) = 0 otherwise. Therefore any joint probability terms with fields which don't contain  will vanish from Equation <ref>.


We can simplify Equation (<ref>) for some specific circumstances:



  *  is not located on any fields. (Θ_i  |) = 0   ∀   i. All terms in the expansion are 0
	
    (⋃_i=1^M _i) = 0


  *  is located on only one patch, A. 
(Θ_i  |) = 0   ∀   i≠ A, (Θ_A |) = 1
	
    (⋃_i=1^M _i) = (_A |Θ_A, 𝐯)


  *  is on the intersection between two fields denoted by A and B:
	
    (⋃_i=1^M _i) =    (_A |Θ_A, 𝐯) + (_B |Θ_B, 𝐯)
       - (_A |Θ_A, 𝐯) ×(_B |Θ_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, λ

    (v_i | n=1, λ) =λ(v_i)/∫dv λ(v)

Expanding this to observations of a population of N objects

    (v_1, v_2...v_N | n=N, λ) = ∏_i=1^N λ(v_i)/∫dv λ(v)

The likelihood of the data is then given by

    (v_1, v_2...v_N, n=N    |λ) 
    
    =   (v_1, v_2...v_N | n=N, λ)(n=N |λ) 
    
    =   ∏_i=1^N λ(v_i)/(∫dv λ(v))^N(∫dv λ(v))^N exp(-∫dv λ(v))/N!
    
    =   ∏_i=1^N λ(v_i) exp(-∫dv λ(v))/N!

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

    ({v}|λ) ∝∏_i=1^N λ(v_i) exp(-∫dv λ(v))

And taking the log likelihood

    ln({v}|λ) ∝ -∫dv λ(v) + ∑_i=1^N lnλ(v_i)





§ MODEL STRUCTURE.



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




In Figure <ref> 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.








