A Flexible Method for Estimating Luminosity Functions via Kernel Density Estimation

Let denote a d-dimensional random vector with density f(x) defined on R^d, and let {} be a sequence of independent identically distributed random variables drawn from f(x). The general d-dimensional multivariate KDE of f is given by

where K( · ) is a multivariate kernel function, and H is the d × d bandwidth matrix or smoothing matrix, which is symmetric and positive definite. In this paper, we focus on the two-dimensional (2D) case.

KDE can be viewed as a weighted sum of density “bumps” that are centered at each data point (Gramacki 2018). The shape of the bumps is determined by the kernel function K( · ), while the width and orientation of the bumps are determined by the bandwidth matrix . It is widely accepted that the choice of kernel function is a secondary matter in comparison with the selection of the bandwidth (e.g., Wand & Jones 1995; Botev et al. 2010; Chen 2017; Gramacki 2018). In most cases, the kernel has the form of a standard multivariate normal density given by

which is what we shall use hereafter. Choosing the form of H depends on the complexity of the underlying density (Sain 2002). Depending on the complexity, H can be isotropic, diagonal, or a full matrix, such as

The top panel of Figure 1 provides a toy example of using the KDE (assuming H is diagonal) in 10 hypothetical observations. Ellipses delineate the kernel spread one bandwidth (the semiaxes of each ellipse are h₁ and h₂) away from each hypothetical point. If H is isotropic, these ellipses will degenerate into circles; else, if H is a full matrix, these ellipses will also have an orientation determined by h₁₂. Sain (2002) argued that for most densities, in particular unimodal ones, allowing different amounts of smoothing for each dimension (diagonal bandwidth matrix) is adequate. Indeed, for the problem in this work, we find that the choice of diagonal bandwidth matrix is typically sufficient. If H is diagonal, Equation (1) can be simplified to

Zoom In Zoom Out Reset image size

Equation (1) assumes that the bandwidth H is constant for every individual kernel. This may lead to poor estimator performance for heterogeneous point patterns (Davies & Baddeley 2018). A popular solution to this problem is using a variable bandwidth or adaptive kernel estimator, which allows H to vary depending on the local density of the input data points (Breiman et al. 1977); a relatively small bandwidth is needed where observations are densely distributed, and a large bandwidth is required where observations are sparsely distributed (e.g., Hu et al. 2012). As a visual comparison to fixed-bandwidth estimation, the bottom panel of Figure 1 shows the kernel spread one bandwidth away from each hypothetical point for the adaptive estimator.

For the diagonal bandwidth matrix case, the bivariate adaptive kernel estimator (see Sain 2002) is given by

Abramson (1982) suggested that the bandwidths are inversely proportional to the square root of the target density itself. Obviously, the target density is unknown. According to Silverman (1986), the practical use of Abramson's (1982) approach is implemented as

where β ≡ 1/2, and is a pilot estimate of the unknown density calculated via Equation (3) with a fixed pilot bandwidth matrix h₁₀ and h₂₀ are the overall smoothing parameters for the variable bandwidths referred to as the global bandwidths (Davies et al. 2018).

The differential LF of a sample of astrophysical objects is defined as the number of objects per unit comoving volume per unit luminosity interval,

where z denotes redshift and denotes the luminosity. Due the typically large span of the luminosities, it is often defined in terms of (e.g., Cara & Lister 2008),

where denotes the logarithm of luminosity.

In an actual survey, only a very limited number of objects in the universe can be observed, as our survey only covers a fraction of the sky and is subject to a selection function. Thus, the estimation of LFs is inevitably based on a truncated sample of objects. The integral of the LF over the survey region W should be approximately the sample size n, for sufficiently large n, i.e.,

where Ω is the solid angle subtended by the survey, and dV/dz is the differential comoving volume per unit solid angle (Hogg 1999). The LF is related to the probability distribution of (z, L) as

Once we have an estimate of p(z, L) by the KDE, we can easily convert this to an estimate of ϕ(L, z) using Equation (9). It needs to be emphasized that the domain of p(z, L) is the survey region W, thus the domain of the LF estimated by Equation (9) is also limited to W.

Choosing the bandwidth parameter is the most crucial issue in using KDE to estimate p(z, L), as the accuracy of KDE depends very strongly on the bandwidth matrix (e.g., Gramacki 2018). Basically, all common approaches to bandwidth selection, which measures the closeness of to its target density f, are based on some appropriate error criteria. One of the most well-known error measurements is the mean integrated square error (MISE). The optimal bandwidth can be obtained by minimizing the Asymptotic MISE (AMISE), the dominant quantity in the MISE (Chen 2017). The three major types of bandwidth selectors include the rule of thumb, cross-validation (CV) methods, and plug-in method. In this work, we focus on the CV method, and our selector is based on the likelihood cross-validation (LCV) criterion.

3.2.1. Likelihood Cross-validation

The LCV thinks about the estimated density itself as a likelihood function, and the LCV objective function is given by

where is the leave-one-out estimator

If H is diagonal, Equation (11) can be simplified to

The LCV bandwidth matrix is the maximizer of the LCV(H) objective function (e.g., Davies et al. 2018)

and it has to be maximized numerically. In view of the difficult numerical procedure as the dimension of data increases, Zhang et al. (2006) improved the LCV method from a Bayesian perspective. They treated nonzero components of H as parameters, whose posterior density can be obtained using the MCMC technique based on the LCV criterion. Given H, the logarithmic likelihood function of {} is

However, as pointed out by Zhang et al. (2006), when using the likelihood given by Equation (15) to perform the Bayesian inference, sufficient prior information on the components of H is required. Using this likelihood function in our problem does not lead to convergence for the MCMC algorithm; we therefore consider a new likelihood function.

3.2.2. A New LCV Method

The new candidate likelihood function is derived from Marshall et al. (1983). In their analysis, the probability distribution of the likelihood for the observables is described by Poisson probabilities. The L−z space is divided into small intervals of size dLdz. In each element, the expected number of objects at z, L is . The likelihood function, , is defined as the product of the probabilities of observing one object at each () element and the probabilities of observing zero objects in all other (z_j, L_j) elements in the accessible regions of the z−L plane (Marshall et al. 1983). Using Poisson probabilities, one has

Transforming to the standard expression and dropping terms that are not model dependent, we obtain

The likelihood function given above naturally incorporates the selection function of the sample (by considering the integral interval W), and it has proved popular (Johnston 2011) and has been widely applied to the parametric estimation of LFs (e.g., Ajello et al. 2012; Yuan et al. 2016; Kulkarni et al. 2019).

Inserting Equation (9) into Equation (16) and dropping terms that are not model dependent, we obtain

Equation (17) enables one to transform the LF measurement into a problem of probability density estimation. can be replaced with its KDE form given by Equation (3), and p(z_i, L_i) is replaced with the leave-one-out estimator constructed by Equation (12). Then, we have

Note that in the first term should not be confused with . The term “leave one out” means that the observations used to calculate are independent of (z_i, L_i). This is where the name “cross-validation” comes from: using one subset of data to make analysis on another subset (see Gramacki 2018).

Using Equation (18), we can obtain the optimal bandwidth by minimizing S,

We find that the above bandwidth selector performs better than the LCV selector given by Equation (13). We note that the numerical optimization result of the LCV selector depends on the initial values, suggesting it may have more than one local maximum in the objective function. In this respect, our new selector has better stability.

Equation (18) can also be used to perform Bayesian inference by combining it with prior information on the bandwidths h₁ and h₂. We find that using uniform priors for h₁ and h₂ are sufficient to employ the MCMC algorithm. Therefore, our new likelihood function is a better choice than the one given by Equation (15), which requires more critical prior information.

Recall that the domain of in Equation (9) is the survey region W. This means that the estimate of using KDE is based on bounded data, while certain difficulties can arise at the boundaries and near them, known as boundary effects or boundary bias (e.g., Müller & Stadtmüller 1999).

In astronomy, many surveys are flux limited. When the observation points are plotted on the L−z plane, there is a clear truncation bound (see Figure 2(a)) defined by f_lim(z):

where d_L(z) is the luminosity distance, F_lim is the survey flux limit, and K(z) represents the K correction. One has for a power-law emission spectrum of index α (e.g., Singal et al. 2011). To give reliable results, any LF estimator must treat the truncation boundary, or flux limit of the survey, properly. From Equation (20), we can see that only when and α is constant is the truncation bound defined by f_lim(z) a curve. Otherwise, if the K correction K(z) takes a more complex form, the truncation bound will be a 2D region but not a simple curve. In this section, we only consider the simple case for the truncation boundary—a simple curve.

Zoom In Zoom Out Reset image size

The mathematical description for the above is that, suppose we observe n points in 2D space, . These points are always from a bounded subset W of the plane, ; W is referred to as the study window (e.g., Davies et al. 2018) or survey region. Direct use of the KDE formula in W may lead to the problem of underestimating the density at boundaries. The reason for the boundary problem is that the kernels from data points near the boundary lose their full probability weight and points that lie just outside the boundary have no opportunity to contribute to the final density estimation on W (e.g., Davies et al. 2018). The boundary effect has long been recognized (Gasser & Müller 1979), and it is an import research topic in KDE study (e.g., Marron & Ruppert 1994; Hall & Park 2002; Marshall & Hazelton 2010; Maleca & Schienle 2014). The common techniques for reducing boundary effects include reflection of data (e.g., Jones 1993), transformation of data (e.g., Marron & Ruppert 1994), and boundary kernel estimators (e.g., Hall & Park 2002). For the LF estimate problems, we find that the transformation method and a method of combining transformation and reflection (Karunamuni & Alberts 2005) are able to reduce boundary effects.

3.3.1. The Transformation Method

The basic idea is as follows: take a one-to-one and continuous function g, and use the regular kernel estimator with the transformed data set . For the LF estimate, we use the following transformation:

where f_lim(z) is given by Equation (20), and δ₁ and δ₂ are transformation parameters. Figure 2(b) provides a toy example showing how a typical flux-limited data set will look like after transformation by Equation (21).

The Jacobian matrix for the above transformation is

and . After the transformation, the density of (x, y), denoted as , can be estimated from Equation (3), and its leave-one-out estimator, , is constructed from Equation (12). Then, by transforming back to the density of the original data set, we have

Inserting Equation (22) into (18), we obtain the negative logarithmic likelihood function S. Following Liu et al. (2011), we estimate the optimal bandwidths and transformation parameters δ₁ and δ₂ simultaneously. They can be estimated by numerically minimizing the object function S, or by combining with uniform priors for h₁, h₂, δ₁, and δ₂, one can also employ the MCMC algorithm to sample the bandwidth and transformation parameters simultaneously. The MCMC algorithm used in this work is “CosmoMC,” a public Fortran code of Lewis & Bridle (2002). Once we know the probability density , the KDE of the LF is easily obtained by

3.3.2. The Transformation–Reflection Method

Introduced by Karunamuni & Alberts (2005), the transformation–reflection method was originally used for univariate data. We extend the method to the bivariate case, and extending it to the trivariate case is also possible (see Section 6.2). First, we transform the original data by

where the meaning of f_lim and δ₁ is similar to those in Equation (21). The determinant of the Jacobian matrix for the above transformation is . Second, we add the missing “probability mass” (e.g., Gramacki 2018) represented by the data set . To illustrate the above ideas, Figure 1(c) shows how a typical flux-limited data set will look after transformation and reflection. The KDE of the density of (x, y) is

The corresponding leave-one-out estimator is

Then, by transforming back to the density of the original data set, we have

Inserting Equations (27) into (18), we can obtain the likelihood function for the transformation–reflection estimator. The next steps for estimating the bandwidth and transformation parameters are similar to those described in Section 3.3.1. Finally, we can obtain the LF estimated by the transformation–reflection approach,

3.3.3. The Transformation–Reflection Adaptive KDE Approach

Figure 3 shows the RLFs for the simulated sample, estimated by the transformation KDE (denoted as ) and the transformation–reflection KDE (denoted as ) methods at several redshifts. The optimal bandwidth and transformation parameters in our KDE methods are summarized in Table 1. Their posterior probability distributions and 2D confidence contours are given in Figure 5. The regular unimodal feature of each distribution suggests that all the parameters are very well constrained (e.g., Yan et al. 2013).

Zoom In Zoom Out Reset image size

Table 1.  Parameters of KDE for the Simulated Sample

Method	Parameters

Note. is implemented on the basis of , and , , and inherit the parameter values of . Parameter errors correspond to the 68% confidence level. Parameters without an error estimate were kept fixed during the fitting stage.

Download table as:  ASCIITypeset image

For comparison, the result measured by the traditional binning method (denoted as ) of Page & Carrera (2000) is also shown. In the literature, bins are commonly chosen somewhat arbitrarily. Yuan & Wang (2013) argued that this may lead to significant bias for the LF estimates close to the flux limit (see panel (a) of Figure 4). In this work, we use the simple rule of thumb suggested by Yuan & Wang (2013) to divide bins. In panel (b) of Figure 4, we illustrate this scheme for dividing bins.

Zoom In Zoom Out Reset image size

Figure 3 shows that both and are generally superior to the binning method, especially for the high-redshift (z ≥ 3) LF estimation. The KDE approaches (especially ) produce relatively smooth LFs, while the binning method produces discontinuous estimates, being prone to have artificial bulges and hollows.

In addition, the KDE approaches have the advantage that their measurement can be extrapolated appropriately beyond the observational limits. In Figure 3, the vertical dashed lines mark the higher luminosity limits of the simulated survey. Note that even beyond the limits, the extrapolated LFs can be a generally acceptable approximation to the true LF.

and have their own advantages and disadvantages. makes smother estimates, but it is prone to producing a slight negative bias at the faint end of the LFs. This is a typical boundary bias, suggesting that a simple transformation method cannot fully solve the boundary problem. performs well at the faint end, but not ideally at the bright end of the LFs. This is expected to be improved by using the adaptive KDE, where the data density is used.

Zoom In Zoom Out Reset image size

Figure 6 shows the RLF and its uncertainty estimated by the transformation–reflection adaptive KDE approach (denoted as ) at several redshifts for the simulated sample. Given posterior probability distributions for the bandwidth and transformation parameters in our KDE methods (see Figure 5), the uncertainty regions of the estimated RLF are plotted by “fgivenx,” a public Python package of Handley (2018). We find that performs well for all six redshifts, and it achieves an excellent approximation to the true LF. Comparing with Figure 3, we notice that can overcome the shortcomings of and .

Zoom In Zoom Out Reset image size

In order to rule out the possibility that our adaptive KDE estimator only gives a good result by chance, we simulate 200 flux-limited samples with their flux limits randomly drawn between 10^−2.5 and 10^−0.5 Jy. All simulated samples share the same input LF (the model A RLF of Yuan et al. 2017). By adjusting the simulated solid angle, we can control the size (n) of each sample. We let the sample size be linearly proportional to the logarithm of its flux limit. Finally, for the 200 simulated samples, their sizes are randomly distributed between 2000 and 40,000 sources.

Figures 7 and 8 shows the RLFs estimated by and , respectively, based on 200 simulated samples. Obviously, the estimator shows better stability and reliability in performance than the binning method. In most cases, we find that the LFs estimated by look like an irregular sawtooth, randomly leaping up and sloping down. This may mislead the observer to use the wrong parametric form to model the LF. Our estimator does not have such a drawback. In fact, this is why the KDE method is very popular in modern statistics: it can produce a smoother estimation that converges to the true density faster (Wasserman 2006).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that in Figure 7 the adaptive KDE method is slightly asymmetrically biased to values above the true LF compared to below the true LF (also see Figure 6), especially for the tails where LFs are extrapolated. This is because the functions of adaptive kernels given in Equation (5) only have one adjustable parameter, β, and cannot let the estimated LF go too steep too soon in the extremely low number density regime. The above effect does not matter as it is only visible for the extrapolated LFs.

To quantify the performance of each LF estimator, we define a statistic d_LF which measures the discrepancy of the estimated LF from the true LF ϕ. If a large number of random vectors, denoted by {}, can be drawn from , can be estimated (e.g., Zhang et al. 2006) as

For the estimator, estimates are given only at some discontinuous points, i.e., the centers of each bin. The statistic d_LF is estimated by

where N_bin is the number of bins, and locates the center of jth bin. The edges of the redshift bins are 0.0, 0.2, 0.5, 0.8, 1.2, 1.8, 2.2, 2.7, 3.3, 3.8, 4.2, 4.7, 5.3, and 6.0. The logarithmic luminosity bins are in increments of 0.3.

In Table 2, we show the statistic d_LF of the different LF estimators calculated for our simulated sample. Also, the mean values for 200 simulated samples are given. The statistic d_LF of all our KDE estimators, especially , is significantly better than that of . d_LF can be understood as the typical error of an LF estimator. In this sense, the estimator improves the accuracy by nearly an order of magnitude compared to . Figure 9 shows the distributions of d_LF of different LF estimators for 200 simulated samples. The smaller dispersion of distributions for our KDE methods suggests that they all have significantly better stability than the estimator, as expected.

Zoom In Zoom Out Reset image size

Table 2.  The Statistic d_LF

dLF	0.0944	0.0239	0.0193	0.0157
	0.1387	0.0311	0.0305	0.0227

Note. d_LF is calculated for the simulated sample described in Section 4. is the mean value for 200 simulated samples.

Download table as:  ASCIITypeset image

Hitherto, the classical binned estimator () is still the most popular nonparametric method. The original version of is the famous 1/V_max estimator (Schmidt 1968). It was originally a variation of the V/V_max test (Schmidt 1968), while V/V_max is essentially a completeness estimator (see Johnston 2011). This makes it impossible for the 1/V_max method to be a mathematically rigorous density estimator. The estimator cannot accurately calculate the surveyed regions for objects close to the flux limit of their parent sample and thus produces a significant systematic error (see Page & Carrera 2000; Yuan & Wang 2013). To improve the 1/V_max estimator, Page & Carrera (2000) proposed a new method that rests on a stronger mathematical foundation. The new method is superior to the 1/V_max estimator (e.g., Yuan & Wang 2013). For the LF at the center of a bin with luminosity interval (L_min, L_max) and redshift interval (z_min, z_max), the new estimator gave

where N is the number of sources detected within the bin. The double integral in Equation (33) naturally considers the truncation boundary as discussed in Section 3.3. In this text, we do not strictly distinguish among the Page & Carrera (2000) method, 1/V_max, and 1/V_a (the generalized version of 1/V_max used for multiple samples; Avni & Bahcall 1980). They are collectively referred to as .

is widely acknowledged for its simplicity and ease of implementation. However, it has some severe drawbacks. First, its estimate depends on the division of bins but currently there are no effective rules to guide the binning. Second, it produces discontinuous estimates, and the discontinuities of the estimate are not due to the underlying LF but are only an artifact of the chosen bin locations. These discontinuities make it very difficult to grasp the structure of the data. Third, the binning of data undoubtedly leads to information loss, and potential biases can be caused by evolution within the bins. Fourth, is a bivariate estimator, and it cannot be properly extended to the trivariate situation. Because the number of bins grows exponentially with the number of dimensions, in higher dimensions one would require many more data points or else most of the bins would be empty (the so-called curse of dimensionality; Gramacki 2018). A trivariate LF estimator is necessary in the situation when additional quantities (such as the photon index; Ajello et al. 2012) besides the redshift and luminosity are incorporated into the LF analysis to tackle complex K corrections (see Yuan et al. 2016). In this case, there is the risk of using a low-dimensional tool to deal with higher dimensional problems. These drawbacks prevent from being a precise estimator. Usually, the estimates of are used to provide guidance or calibration for modeling the LF. If the result itself is not accurate, the parametric description cannot be reliable.

From a mathematical perspective, the estimator is a kind of 2D histogram. All the drawbacks of are inherited from the histogram. In the mathematical community, the histogram is an outdated density estimator, used only for rapid visualization of results in one or two dimensions (e.g., Gramacki 2018). To conquer the shortcomings of histograms, mathematicians have developed some other nonparametric estimators, such as smoothing histograms, Parzen windows, k-nearest neighbors, KDE, etc. Among these, KDE is the most popular density estimator. This is why we develop the new LF estimator within a KDE framework. Our KDE methods can overcome all the drawbacks of and significantly improve the accuracy.

The Lynden-Bell (1971) method, and its variants (e.g., Efron & Petrosian 1992; Caditz & Petrosian 1993), is another important estimator. It overcomes some of the shortcomings in . Recently, Singal et al. (2014) extended C⁻ to the trivariate scenario and presented the redshift evolutions and distributions of the gamma-ray luminosity and photon spectral index of flat-spectrum radio quasars. This suggested that C⁻ has good extensibility in dimension. However, the direct product of C⁻ is the cumulative distribution function of the LF but not the LF itself, and this is inconvenient. A more serious disadvantage is that typically assumes that luminosity and redshift are statistically independent. In actual application, one inevitably has to introduce a parameterized model that removes the correlation of luminosities and redshifts (e.g., Singal et al. 2014).

Motivated by the potential hazards and pitfalls existing in the above traditional estimators, approaches based on more innovative statistics have been developed. One typical representative is the semiparametric approach of Schafer (2007). Schafer decomposed the bivariate density into

where f(z) and g(z) are determined nonparametrically, while h(z, M, θ) has a parametric form. From Figure 11, we find that our adaptive KDE method achieves estimates comparable to those of Schafer (2007) while making far fewer assumptions about the shape of the LF. The Bayesian approach of Kelly et al. (2008) is another estimator based on innovative statistics. This method is similar to ours in some aspects, e.g., both are within a Bayesian framework and using MCMC algorithms for estimating the parameters; both model the LF as a mixture of Gaussian functions. The difference is that we arrange a relatively unified Gaussian function at each data point, and the final LF is the sum of all the Gaussian functions, while Kelly et al. (2008) use far fewer (typically ∼3–6) Gaussian functions, and each of them has six free parameters for adjusting the shape and location. Therefore, our method has much fewer (typically ∼3–6) parameters than theirs (typically ∼22–40). In addition, their method requires much more critical prior information for parameters to aid the convergence of the MCMC. Generally, they assume that the LF is unimodal, thus their prior distributions are constructed to place more probability on situations where the Gaussian functions are close together (see Figure 3 of Kelly et al. 2008). Our KDE method does not impose any preliminary assumptions on the shape of LFs, and simply use uniform (so-called “uninformative”) priors for the parameters to run the MCMC. Therefore, our method is more flexible.

Our KDE methods can be easily extended to the trivariate LF situation. Below we give an example for employing the estimator in this case. Suppose we observe n objects within a survey region W, and . are the single power-law spectral index, redshift, and luminosity of the ith object, respectively. Transform the observed data by

where defines the truncation boundary of the sample and is given by

where d_L(z) is the luminosity distance, is the survey flux limit, and represents the K correction, where for a power-law emission spectrum of index α. If the emission spectrum is not a power law, would have a different functional form. Thus, we reserve the possibility that α can also represent other quantities to determine the K correction. In Appendix C, we give the procedure of using the transformation–reflection KDE for trivariate LFs. The adaptive KDE formulas can be easily derived.

In mathematics, KDE is a very effective smoothing technique with many practical applications, and analysis of KDE methods is an ongoing research task. New research on bandwidth selection, reducing the boundary bias, multivariate KDE, and fast KDE algorithm for big data, etc. are continuously emerging (e.g., Cheng et al. 2019; Igarashia & Kakizawa 2019). The fast KDE algorithm is especially relevant to astronomical surveys in the near future. From Equation (1), we note that the computational complexity of KDE is of , and it will be computationally expensive for large data sets and higher dimensions. In recent years, techniques such as using the fast Fourier transform have been proposed to accelerate the KDE computations (e.g., Gramacki & Gramacki 2017, 2017; Davies & Baddeley 2018). These ongoing developments in KDE theory mean that our LF estimator could be seen as a starting method with room to upgrade and enable it to flexibly deal with various LF-estimating problems in future surveys.