Tools for Conditional Density Estimation

Conditional Density Estimation Tools in Python and R
with Applications to Photometric Redshifts and Likelihood-Free Cosmological Inference


It is well known in astronomy that propagating non-Gaussian prediction uncertainty in photometric redshift estimates is key to reducing bias in downstream cosmological analyses. Similarly, likelihood-free inference approaches, which are beginning to emerge as a tool for cosmological analysis, require a characterization of the full uncertainty landscape of the parameters of interest given observed data. However, most machine learning (ML) or training-based methods with open-source software target point prediction or classification, and hence fall short in quantifying uncertainty in complex regression and parameter inference settings such as the applications mentioned above. As an alternative to methods that focus on predicting the response (or parameters) from features , we provide nonparametric conditional density estimation (CDE) tools for approximating and validating the entire probability density function (PDF) of given (i.e., conditional on) . This density approach offers a more nuanced accounting of uncertainty in situations with, e.g., nonstandard error distributions and multimodal or heteroskedastic response variables that are often present in astronomical data sets. As there is no one-size-fits-all CDE method, and the ultimate choice of model depends on the application and the training sample size, the goal of this work is to provide a comprehensive range of statistical tools and open-source software for nonparametric CDE and method assessment which can accommodate different types of settings — involving, e.g., mixed-type input from multiple sources, functional data, and images — and which in addition can easily be fit to the problem at hand. Specifically, we introduce four CDE software packages in Python and R based on ML prediction methods adapted and optimized for CDE: NNKCDE, RFCDE, FlexCode, and DeepCDE. Furthermore, we present the cdetools package with evaluation metrics. This package includes functions for computing a CDE loss function for tuning and assessing the quality of individual PDFs, together with diagnostic functions that probe the population-level performance of the PDFs. We provide sample code in Python and R as well as examples of applications to photometric redshift estimation and likelihood-free cosmological inference via CDE.

Nonparametric statistics, Statistical software, Statistical computing, methods: data analysis, galaxies: distances and redshifts, cosmology: cosmological parameters

Niccolò Dalmasso,

1 Introduction

Machine learning (ML) has seen a surge in popularity in almost all fields of astronomy that involve massive amounts of complex data (Ball & Brunner, 2010; Way et al., 2012; Ivezić et al., 2014; Ntampaka et al., 2019). Most ML methods target regression and classification, whose primary goal is to return a point estimate of an unknown response variable given observed features , often falling short in quantifying nontrivial uncertainty in . For instance, returning a point estimate for a supernova’s type given a supernova’s light curve , or for a galaxy mass given its light spectrum , fails to capture degeneracies in the mapping from to . Neglecting to propagate these uncertainties through down-stream analyses may lead to imprecise or even inaccurate inferences of physical parameters.

Consider the following two examples of problems where uncertainty quantification can be impactful:

  • Photometric redshift estimation. In photometric redshift (photo-) estimation, one attempts to constrain the cosmological redshift () of a galaxy after observing the shifted spectrum using a handful of broadband filters, and sometimes additional variables such as morphology and environment. In a prediction setting, the response could be the galaxy’s true (i.e. spectroscopically observed) redshift but could also include other galaxy properties () such as the galaxy’s mass or age; the features would represent the collection of directly observable inputs such as photometric magnitudes and colors used to predict or more generally a multivariate response . The use of photo- posterior estimates — that is, estimates of the probability density functions (PDFs) for individual galaxies — is crucial for cosmological inference from photometric galaxy surveys, as of currently cataloged galaxies are observed solely via photometry, a percentage that will only grow in the coming decade as the Large Synoptic Survey Telescope (LSST) begins gathering data (Ivezić et al., 2019). In addition, widely different redshifts can be consistent with the observed colors of a galaxy; photo- posterior estimates can capture such degeneracy or multimodality in the distribution whereas point estimates cannot. Though the benefits of using posteriors over have yet to be fully exploited (Viironen et al., 2018), it is thoroughly established that one can improve down-stream cosmological analysis by properly propagating photo- estimate uncertainties via probability density functions (PDFs) rather than just using simple point predictions of (Mandelbaum et al., 2008; Wittman, 2009; Sheldon et al., 2012; Carrasco Kind & Brunner, 2013; Graff et al., 2014; Schmidt et al., 2019).

  • Forward-modeled observables in cosmology. Some cosmological probes, such as the type Ia supernova (SN Ia) distance-redshift relationship, have observables that are straightforward to simulate in spite of an intractable likelihood. Likelihood-Free Inference (LFI) methods allow for parameter inference in settings where the likelihood function, which relates the observable data to the parameters of interest , is too complex to work with directly, but one is able to simulate from a stochastic forward model at fixed parameter settings . The most common approach to LFI or simulation-based parameter inference is Approximate Bayesian Computation (ABC), whose many variations (see Beaumont et al. 2002 and Sisson et al. 2018 for a review) repeatedly draw from the simulator and compare the output with observed data in real time to arrive at a set of plausible parameter values consistent with . With computationally intensive simulations, however, a classical ABC approach may not be practical. An alternative approach to ABC rejection sampling is to use faster training-based methods to assess the uncertainty about for any first, and then consider the specific case .

From a statistical perspective, the right tool for quantifying the uncertainty about once is observed is the conditional density . In a prediction context such as for photo- problems, where heteroskedastic errors or multimodal responses may occur, conditional density estimation (CDE) of the density for the redshift of individual galaxies given photometric data provides a more nuanced accounting of uncertainty than a point estimate or prediction interval alone. CDE can also be used in LFI where, in our notation, the parameters of interest take the role of the “response” . In such settings, one can apply training-based approaches to forward-simulated data to estimate the posterior probability distribution of cosmological parameters given observed data , and from these posteriors then derive, e.g. posterior credible intervals of . Section 4.2 shows an example of cosmological parameter inference using mock weak lensing data. Here we follow Izbicki et al. (2014, 2019) to combine ABC and CDE by directly applying CDE techniques (Section 2) and loss functions (Section 3.1) to simulated data . Similar works include performing LFI via CDE using Gaussian copulas (Li et al., 2017; Chen & Gutmann, 2019) and random forests (Marin et al., 2016). Other examples include neural density estimation in LFI via mixture density networks and masked autoregressive flows (Papamakarios & Murray, 2016; Lueckmann et al., 2017, 2019; Papamakarios et al., 2019; Greenberg et al., 2019; Alsing et al., 2018, 2019).

Data in astronomy present a challenge to estimating conditional densities, due to both the complexity of the underlying physical processes and the complicated observational noise. Precision cosmology, for example, requires combining data from different scientific probes, each affected by unique sources of systematic uncertainty, to produce samples from complicated joint likelihood functions with nontrivial covariances in a high-dimensional parameter space (Krause et al., 2017; Joudaki et al., 2017; Aghanim et al., 2018; van Uitert et al., 2018; Abbott et al., 2019). In such situations, CDE methods that target a variety of settings and non-standard data (images, correlation functions, mixed data types) become especially valuable. However, for any given data type, there is no one-size-fits-all CDE method. For example, deep neural networks often perform well in settings with large amounts of representative training data but in applications with smaller training samples one may need a different tool. There is also additional value in models that are interpretable and easy to fit to the data at hand.

The goal of this paper is to provide statistical tools and open-source software for nonparametric CDE and method assessment appropriate for challenging data in a variety of inference scenarios.

To our knowledge, existing CDE software either targets discrete (e.g., probabilistic classifiers or ordinal classification (Frank & Hall, 2001)) or uses kernel density estimation (KDE) across all data points to provide an estimate for a continuous (Hyndman et al., 2018). What distinguishes our methodology work from others is that we are able to adapt virtually any training-based prediction method to the problem of estimating full probability distributions. By leveraging existing ML methods, we are hence able to provide uncertainty prediction methods and software for more general and complex data settings than ’the computational tools currently available in the literature. In this paper, we showcase and provide Python and R code for four flexible CDE methods: NNKCDE, RFCDE/fRFCDE, FlexCode and DeepCDE.1 Each CDE approach has particular usage for different settings of response dimensionality, feature types, and computational requirements. Table 1 provides a high-level overview in the top panel, and lists some properties of each method in the bottom panel.

A highlight of our software is that it makes uncertainty quantification straightforward for users of standard open-source machine learning Python packages. As NNKCDE, RFCDE/fRFCDE, FlexCode share the sklearn (Pedregosa et al., 2011) API (with fit and predict methods), our methods are usable within the sklearn ecosystem for, e.g., cross-validation and model selection. In addition, FlexCode is essentially a plug-in method where the user can utilize any sklearn-compatible regression function. DeepCDE has implementations for both Tensorflow (Abadi et al., 2015) and Pytorch (Paszke et al., 2019), two of the most widely used deep learning frameworks.

Method Name Summary Hyper-parameters Details
NNKCDE Computes a KDE estimate of
Nearest Neighbor multivariate using the nearest Number of neighbors Section 2.1
Kernel CDE neighbors of the evaluation point Kernel bandwidth (Code: A.1)
in feature space.
RFCDE Random forest that partitions
the feature space using a CDE loss. Random forest hyperparams. Section 2.2
Random Forest CDE Constructs a weighted KDE estimate Kernel bandwidth (Code: A.2)
of multivariate with weights
defined by leaves in the forest.
fRFCDE RFCDE version suitable for functional Random forest hyperparams.
functional Random features . Partitions the feature Kernel bandwidth Section 2.2.1
Forest CDE space directly rather than Partition parameter (Code: A.2.1)
representing as a vector.
FlexCode Flexible Conditional Uses basis expansion of univariate Number of expansion coeffs. Section 2.3
Density Estimation to turn CDE into a series of Selected regression method (Code: A.3)
univariate regression problems. hyperparams.
DeepCDE Uses basis expansion of univariate
Deep Neural similar to FlexCode, but learns Number of expansion coeffs. Section 2.4
Networks CDE the expansion coefficients Selected deep neural network (Code: A)
simultaneously using a deep architecture hyperparams.
neural network.



Method Capacity (# Training Pts) Multivariate Response Functional Features Image Features NNKCDE Up to (f)RFCDE Up to FlexCode Up to DeepCDE Up to
Table 1: Top: Naming convention, high-level summary and hyper-parameters of CDE methods, along with references for further details and code examples. Bottom: Comparison of CDE methods in terms of training capacity and compatibility with multivariate response and different types of features, with capacities estimated based on input with around 100 features and a standard i5/i7/quad-core processor with GB of RAM. Note that less complex methods (such as NNKCDE) tend to be easier to use, easier to interpret, and often perform better in settings with smaller training sets, whereas more complex methods (such as DeepCDE) perform better in settings with larger (representative) training sets.

In addition to the CDE methods above, we provide the package cdetools, which can be used for tuning and assessing the performance of CDE models on held-out validation data. CDE method assessment is challenging per se because we never observe the true conditional probability density, merely samples (observations) from it. Furthermore, whereas loss functions such as the root-mean-squared error (RMSE) are typically used in regression problems, they are not appropriate for the task of uncertainty quantification of estimated probability densities. The cdetools package provides two types of functions for method assessment. First, it provides functions for computing a so-called CDE loss function (defined by Equation 4 in Section 3.1) for tuning and assessing the quality of individual PDFs. Second, it provides diagnostic functions that probe the population-level performance of the PDFs. More specifically, we have included functions for computing the Probability Integral Transform (PIT) and Highest Posterior Density (HPD); these metrics check how well the final density estimates on average fit the data in the tail and highest-density regions, respectively (see Section 3.2, and Figure 2 for a visual sketch).

Organization of the paper. In Section 2, we introduce tools for conditional density estimation (NNKCDE, RFCDE/fRFCDE, FlexCode, DeepCDE). In Section 3, we describe tools for model selection and diagnostics. Then, in Section 4, we illustrate our CDE and method assessment tools for three different applications: photo- estimation, likelihood-free cosmological inference and spec- estimation. Python and R code usage examples can be found in Appendix A.

Notation. We denote the true (unknown) conditional density by , whereas an estimate of the density is denoted by . We represent the CDE loss function that measures the discrepancy between the conditional density and its estimate by . Typically this loss cannot be computed directly because it depends on unknown quantities; an estimate of the loss is denoted by . As before, we continue to use bold-faced letters to denote vectors.

2 Overview of Conditional Density Estimation Tools

We start by briefly describing the conditional density estimators in Table 1. Unless otherwise stated, we choose the tuning or hyper-parameters by minimizing the CDE empirical loss in Equation 5 using cross-validation.

2.1 Nnkcde

Nearest-Neighbors Kernel CDE (NNKCDE; Izbicki et al. 2017, Freeman et al. 2017) is a simple and easily interpretable CDE method. It computes a kernel density estimate of using the nearest neighbors of the evaluation point . The model has only two tuning parameters: the number of nearest neighbors and the bandwidth of the smoothing kernel in -space. Both tuning parameters are chosen in a principled way by minimizing the CDE loss on validation data.

More specifically, the kernel density estimate of given is defined as


where is a normalized kernel (e.g., a Gaussian function) with bandwidth , is a distance metric, and is the index of the nearest neighbor of . It is essentially a smoother version of the histogram estimator proposed by Cunha et al. (2009) in that it approximates the density with a smooth continuous function rather than by binning.

We provide the NNKCDE2 software in both Python and R (Pospisil & Dalmasso, 2019a), with examples in Section A.1. NNKCDE can accommodate any number of training examples. However, selecting tuning parameters scales quadratically in , the number of nearest neighbors, which prohibits using . Our implementation (Izbicki et al. 2019, Appendix D) is computationally more efficient than standard nearest-neighbor kernel smoothers. For instance, we are able to efficiently evaluate the loss function in Equation 5 on large validation samples by expressing the first integral in terms of convolutions of the kernel function; these have a fast-to-compute closed solution for a Gaussian kernel.

2.2 Rfcde

Random forests (RFs, Breiman 2001) is one of the best off-the-shelf solutions for regression and classification problems. It builds a large collection of decorrelated trees, where each tree is a data-based partition of the feature space. The trees are then averaged to yield a prediction. RFCDE, introduced by Pospisil & Lee (2018), is an extension of random forests to conditional density estimation and multivariate responses. Like NNKCDE, it computes a kernel density estimate of but with nearest neighbor weightings defined by the location of the evaluation point relative to the leaves in the random forest. RFCDE inherits the advantages of random forests in that it can handle mixed-typed data. It also does not require the user to specify distances or similarities between data points, and it has good performance while remaining relatively interpretable.

The main departure from other random forest algorithms is our criterion for feature space partitioning decisions. In regression contexts, the splitting variable and split point are typically chosen so as to minimize the mean-squared-error loss. In classification contexts, the splits are typically chosen so as to minimize a classification error. Existing random forest density estimation methods such as quantile regression forests by Meinshausen (2006) and the TPZ algorithm by Carrasco Kind & Brunner (2013) use the same tree structure as regression and classification random forests, respectively. RFCDE, however, builds trees that minimize the CDE loss (see Equation 5), allowing the forest to adapt to structures in the conditional density; hence overcoming some of the limitations of the usual regression approach for data with heteroskedasticity and multimodality. In addition, RFCDE does not require discretizing the response as in TPZ, thereby providing more accurate results at a lower cost for continuous responses, especially in the case of multivariate continuous responses where binning is problematic. See Pospisil & Lee (2018) for further examples and comparisons.

Another unique feature of RFCDE is that it can handle multivariate responses with joint densities by applying a weighted kernel smoother to . This added feature enables analysis of complex conditional distributions that describe relationships between multiple responses and features, or equivalently between multiple parameters and observables in an LFI setting. Like quantile regression forests, the RFCDE algorithm takes advantage of the fact that random forests can be viewed as a form of adaptive nearest-neighbor method with the aggregated tree structures determining a weighting scheme. This weighting scheme can then be used to estimate the conditional density , as well as the conditional mean and quantiles, as in quantile regression forests (but for CDE-optimized trees). As mentioned above, RFCDE computes the latter density by a weighted kernel density estimate (KDE) in using training points near the evaluation point . These distances are effectively defined by how often a training point belongs to the same leaf node as in the forest (see Equation 1 in Pospisil & Lee 2018 for details).

Despite the increased complexity of our CDE trees, RFCDE still scales to large data sets because of an efficient computation of splits via orthogonal series. Moreover, RFCDE extends the density estimates on new to the multivariate case through the use of multivariate kernel density estimators (Epanechnikov, 1969). In both the univariate and multivariate cases, bandwidth selection can be handled by either plug-in estimators or by tuning using a density estimation loss.

For ease of use in the statistics and astronomy communities, we provide RFCDE3 in both Python and R, which call a common C++ implementation of the core training functions that can easily be wrapped in other languages (Pospisil & Dalmasso, 2019b).

Remark: Estimating a CDE loss is an inherently harder task than calculating the mean squared error (MSE) in regression. As a consequence, RFCDE might not provide meaningful tree splits in settings with a large number of noisy features. In such settings, one may benefit from combining a regular random forest tree structure (optimized for regression) with a weighted kernel density estimator (for the density calculation). See Section 4.3 for an application to functional data along with software implementation.4


In addition the RFCDE package includes fRFCDE, a variant of RFCDE (Pospisil & Lee, 2019), that can accommodate functional features by partitioning in the continuous domain of such features. The spectral energy distribution (SED) of a galaxy is its energy as a function of continuous wavelength of light; hence it can be viewed as functional data. Another example of functional data is the shear correlation function of weak lensing, which measures the mean product of the shear at two points as a function of a continuous distance between those points. Similarly, any function of continuous time is an example of functional data. Treating functional features (like spectra, correlation functions or images) as unordered multivariate vectors on a grid suffers from a curse of dimensionality. As the resolution of the grid becomes finer the dimensionality of the data increases but little additional information is added, due to high correlation between nearby grid points. fRFCDE adapts to this setting by partitioning the domain of each functional feature (or curve) into intervals, and passing the mean values of the function in each interval as input to RFCDE. Feature selection is then effectively done over regions of the domain rather than over single points. More specifically, the partitioning in fRFCDE is governed by the parameter of a Poisson process, with each functional feature entering as a high-dimensional vector . Starting with the first element of the vector, we group the first elements together. We then repeat the procedure sequentially until we have assigned all elements into a group; this effectively partitions the function domain into disjoint intervals . The function mean values or smoothed brightness measurements of each interval are finally treated as new inputs to a standard (vectorial) RFCDE tree. The splitting of the smoothed predictors is done independently for each tree in the forest. Other steps of fRFCDE, such as the computation of variable importance, also proceed as in (vectorial) RFCDE but with the averaged values of a region as inputs. As a result, fRFCDE has the capability of identifying the functional inputs and the regions in the input domain that are important for estimating the response . Figure 1 shows schematically the differences and similarities in construction between standard RFCDE and its fRFCDE variant.

The fRFCDE method is as scalable as standard random forests, accommodating training sets on the order of . As the examples in Section 4.3 show, we can obtain substantial gains with a functional approach, both in terms of statistical performance (that is, CDE loss) and computational time. In addition, the change in code syntax is minimal, as one only needs to pass the Poisson parameter as the flambda argument during the forest initialization. Examples in Python and R are provided in Appendix A.2.1.

Figure 1: A schematic diagram of RFCDE (top row) and fRFCDE (bottom row) applied to a galaxy spectrum from Vanderplas et al. (2012). Top row: RFCDE treats the intensity at each recorded wavelength of the spectrum as a feature or “input” to the random forests algorithm — the blue vertical dashed lines indicate every 100 recorded wavelength. RFCDE then builds an ensemble of CDE trees, where each tree partitions the feature space according to the CDE loss, as illustrated in the top right figure for features and . Bottom row: fRFCDE instead groups nearby measurements together where the group divisions are defined by a Poisson process with parameter (vertical green dashed lines, left figure). The new smoothed features are computed by integrating the intensity over the grouped wavelengths. A forest of CDE trees is then built using the same construction as in RFCDE but with the smoothed features as inputs (bottom right figure).

2.3 FlexCode

Introduced by Izbicki & Lee (2017), FlexCode5 is a CDE method that uses a basis expansion for the univariate response and poses CDE as a series of univariate regression problems. The main advantage of this method is its flexibility as any regression method can be applied towards CDE, enabling us to tailor our choice of regression method to the intrinsic structure and type of data at hand.

More precisely, let be an orthonormal basis like a Fourier or wavelet basis for functions of . The key idea of FlexCode is to express the unknown density as a basis expansion


By the orthogonality property of the basis, the (unknown) expansion coefficients are then just orthogonal projections of onto the basis vectors. We can estimate these coefficients using a training set of data by regressing the transformed response variables on predictors for every basis function (see Izbicki & Lee (2017) Equation 2.2, for details). The number of basis function is chosen by minimizing a CDE loss function on validation data. The estimated density, , may contain small spurious bumps induced by the Fourier approximation and may not integrate to one. We remove such artifacts as described in Izbicki & Lee (2016) by applying a thresholding parameter chosen via cross-validation. FlexCode turns a challenging density estimation problem into a simpler regression problem, where we can choose any regression method that fits the problem at hand.

To provide a concrete example, for high-dimensional (such as galaxy spectra) we can use methods such as sparse or spectral series (“manifold”) regression methods (Tibshirani, 1996; Ravikumar et al., 2009; Lee & Izbicki, 2016); see Section 4.3 for an example with FlexCode-Series. For multi-probe studies with mixed data types, we can use random forests regression (Breiman, 2001). On the other hand, large-scale photometric galaxy surveys such as LSST require methods that can work with data from millions, if not billions, of galaxies. Schmidt et al. (2019) present the results of an initial study of the LSST Dark Energy Science Collaboration (LSST-DESC). Their initial data challenge (“Photo- DC 1”) compares the CDEs of a dozen photo- codes run on simulations of LSST galaxy photometry catalogs in the presence of complete, correct, and representative training data. FlexZBoost, a version of FlexCode based on the scalable gradient boosting regression technique by Chen & Guestrin (2016), was entered into the data challenge because of the method’s ability to scale to massive data. In the DC1 analysis, FlexZBoost was among the the strongest performing codes according to established performance metrics of such PDFs and was one of only two codes to beat the experimental control under a more discriminating metric, the CDE loss.

For massive surveys such as LSST, FlexCode also has another advantage compared to other CDE methods, namely its compact, lossless storage format. Juric et al. (2017) establishes that LSST has allocated 100 floating point numbers to quantify the redshift of each galaxy. As is shown in Schmidt et al. (2019), the myriad methods for deriving photo- PDFs yield radically different results, motivating a desire to store the results of more than one algorithm in the absence of an obvious best choice. For the photo- PDFs of most codes, one may need to seek a clever storage parameterization to meet LSST’s constraints (Carrasco Kind & Brunner, 2014; Malz et al., 2018), but FlexCode is virtually immune to this restriction. Since FlexCode relies on a basis expansion, one only needs to store coefficients per target density for a lossless compression of the estimated PDF with no need for binning. Indeed, for DC1, we can with FlexZBoost reconstruct our estimate at any resolution from estimates of the first 35 coefficients in a Fourier basis expansion. In other words, FlexZBoost enables the creation and storage of high-resolution photo- catalogs for several billion galaxies at no added cost.

Our public implementation of FlexCode— available in both Python (Pospisil et al., 2019) and R (Izbicki & Pospisil, 2019), respectively — cross-validates over regression tuning parameters (such as the number of nearest neighbors in FlexCode-kNN) and computes the FlexCode coefficients in parallel for further time savings. The computational complexity of FlexCode will be the same as parallel individual regressions. In particular, the scalability of FlexCode is determined by the underlying regression method. A scalable method like XGBoost leads to scalable FlexCode fitting. In the Python version of the code, the user can choose between the following regression methods: XGBoost for FlexZBoost  but also nearest neighbors, LASSO (Tibshirani, 1996), and random forests regression (Breiman, 2001). In the R version, the following choices are available: XGBoost, nearest neighbors, LASSO, random forests, Nadaraya-Watson kernel smoothing (Nadaraya, 1964; Watson, 1964), sparse additive models (Ravikumar et al., 2009), and spectral series (“manifold”) regression (Lee & Izbicki, 2016). In both implementations, the user may also use a custom regression method; we illustrate how this can be done with a vignette in both packages. For the Python implementation, the user can add any custom regression method following the sklearn API, i.e., with fit and predict methods. Examples in both languages are presented in Appendix A.3.

2.4 DeepCDE

Recently, neural networks have reemerged as a powerful tool for prediction settings where large amounts of representative training data are available; see LeCun et al. (2015) and Goodfellow et al. (2016) for a full review. Neural networks for CDE, such as Mixture Density Networks (MDNs; Bishop 1994) and variational methods (Tang & Salakhutdinov, 2013; Sohn et al., 2015), usually assume a Gaussian or some other parametric form of the conditional density. MDNs have lately also been used for photometric redshift estimation (D’Isanto & Polsterer, 2018; Pasquet et al., 2019) and for direct estimation of likelihoods and posteriors in cosmological parameter inference (see Alsing et al. 2019 and references within).

DeepCDE6 (Dalmasso & Pospisil, 2019) takes a different, fully nonparametric approach to CDE. It combines the advantages of basis expansions with the flexibility of neural network architectures, allowing for data types like image features and time-series data. DeepCDE is based on the orthogonal series representation in FlexCode, given in Equation 2, but rather than relying on regression methods to estimate the expansion coefficients in Equation 2, DeepCDE computes the coefficients jointly with a neural network that minimizes the CDE loss in Equation 4. Indeed, one can show that for an orthogonal basis, the problem of minimizing this CDE loss is (asymptotically) equivalent to finding the best basis coefficients in FlexCode under mean squared error loss for the individual regressions; see Appendix C for a proof. The value of this result is that DeepCDE with a CDE loss directly connects prediction with uncertainty quantification, implying that one can leverage the state-of-the-art deep architectures for an application at hand toward uncertainty quantification for the same prediction setting.

From a neural network architecture perspective, DeepCDE only adds a linear output layer of coefficients for a series expansion of the density according to


where is an orthogonal basis for functions of . Like FlexCode, we normalize and remove spurious bumps from the final density estimates according to the procedure in Section 2.2 of Izbicki & Lee (2016).

The greatest benefit of DeepCDE is that it can be implemented with both convolutional and recurrent neural network architectures, extending to both image and sequential data. For most deep architectures, adding a linear layer represents a small modification, and a negligible increase in the number of parameters. For instance, with the AlexNet architecture (Krizhevsky et al., 2012), a widely used, relatively shallow convolutional neural network, adding a final layer with 30 coefficients for a cosine basis only adds extra parameters. This represents a increase over the million already existing parameters, and hence a negligible increase in training and prediction time. Moreover, the CDE loss for DeepCDE is especially easy to evaluate; see Appendix C for details.

We provide both TensorFlow and Pytorch implementations of DeepCDE (Dalmasso & Pospisil, 2019). We also include examples that shows how one can easily build DeepCDE on top of AlexNet (Krizhevsky et al., 2012); in this case, for the task of estimating the probability distribution of the correct orientation of a color image . Note that Fischer et al. (2015) use AlexNet for the corresponding regression task of predicting the orientation for a color image but without quantifying the uncertainty in the predictions.

3 How to Assess Method Performance

After fitting CDEs, it is important to assess the quality of our models of uncertainty. For instance, after computing photo- PDF estimates for some galaxies, one may ask whether these estimates accurately quantify the true uncertainty distributions . Similarly, in the LFI task, a key question is whether an estimate of the posterior distribution, of the cosmological parameters is close enough to the true posterior given the observations .

We present two method-assessment tools suitable to different situations, which are complementary and can be performed simultaneously, with public implementations in the cdetools7 package in both Python and R (Dalmasso et al., 2019). First, we describe a CDE loss function that directly provides relative comparisons between conditional density estimators, such as the methods presented in Section 2 or, equivalently, between a set of models (for the same method) with different tuning parameters. Second, we describe visual diagnostic tools, such as Probability Integral Transforms (PIT) and Highest Probability Density (HPD) plots, that can provide insights on the overall goodness-of-fit of a given estimator to observed data.

3.1 CDE loss

Here we briefly review the CDE loss from Izbicki & Lee (2016) for assessing conditional density estimators and discuss it in the context of the cosmology LFI case.

The goal of a loss function is to provide relative comparisons between different estimators, so that it is easy to directly choose the best fitted model among a list of candidates. Given an estimate of , we define the CDE loss by


where is the marginal distribution of the features . This loss is the CDE analog to the standard mean squared error (MSE) in standard regression. The weighting by the marginal distribution of the features emphasizes that errors in the estimation of for unlikely features are less important. The CDE loss cannot be directly evaluated because it depends on the unknown true density . However, one can estimate the loss (up to a constant determined by the true ) by


where represents our validation or test data, i.e., a held-out set not used to construct . In our implementation, the function cde_loss returns the estimated CDE loss as well as an estimate of the standard deviation or the standard error (SE) of the estimated loss.

CDE loss for LFI. In LFI settings, we use a slightly different version of the CDE loss in Eq. 4. Because the goal (in ABC) is to approximate the posterior density at fixed , a natural evaluation metric is the integrated squared error loss


of the conditional density at only. Estimating this loss can however be tricky as only a single instance of data with is available in practice. Hence, Izbicki et al. (2019) approximates Equation 6 by computing the empirical loss in Eq. 5 over a restricted subset of the validation data that only includes the points that fall in an -neighborhood of , where is the tolerance of the ABC rejection algorithm. The detailed analysis of this approach can be found in Izbicki et al. (2019).

3.2 PIT and HPD diagnostics

The CDE loss function is a relative measure of performance that cannot address absolute goodness-of-fit. To quantify overall goodness-of-fit, we examine how well calibrated an ensemble of conditional density estimators are on average, over validation or test data . For ease of notation, we will in this section denote and for a generic by and .

Given a true probability density of a variable conditioned on data , an estimated probability density cannot be well-calibrated unless . Built on the same logic, the probability integral transform (PIT; Polsterer et al. 2016)


assesses the calibration quality of an ensemble of CDEs for scalar representing the cumulative distribution function (CDF) of evaluated at ; this PIT value corresponds to the shaded area in Figure 2, left. A statistically self-consistent population of densities has a uniform distribution of PIT values, and deviations from uniformity indicate inaccuracies of the estimated PDFs. Overly broad CDEs manifest as under-representation of the lowest and highest PIT values, whereas overly narrow CDEs manifest as over-representation of the lowest and highest values.

However, PIT values do not easily generalize to multiple responses. For instance, for a bivariate response , the quantity is not in general uniformly distributed (Genest & Rivest, 2001). An alternative statistic that easily generalizes to multivariate is the highest probability density value (HPD; Izbicki et al. 2017, Appendix A):


The HPD value is based on the definition of the highest density region (HDR, Hyndman 1996) of a random variable ; that is, the subset of the sample space of where all points in the region have a probability above a certain value. The HDR of can be seen as a region estimate of y when is observed. In words, the set is the smallest HDR containing the point , and the HPD value is simply the probability of such a region. Figure 2, right, shows a schematic diagram of the HPD value (green shaded area) and HDR region (highlighted segments on the y-axis) for the estimated density . The HPD value can also be viewed as a measure of how plausible is according to and is directly related to the Bayesian analog of p-values or the e-value (Pereira & Stern, 1999). One can show (Harrison et al., 2015) that even for multivariate , the HPD values for validation data follow a distribution if the CDEs are well-calibrated on the population level. Thus, these values can also be used for assessing the fit of conditional densities in the same way as PIT values. In our implementation the functions pit_coverage(cde_test, y_grid, y_test) and hpd_coverage(cde_test, y_grid, y_test), respectively, calculate PIT and HPD values for CDE estimates cde_test using the grid y_grid with observed values y_test.

The PIT and HPD are not without their limitations, however, as demonstrated in the control case of Schmidt et al. (2019) and Figure 3 of Section 4.1 here. Because the PIT and HPD values can be uniformly distributed even if is not well estimated, they must be used in conjunction with loss functions for method assessment. A popular way of visualizing PIT and HPD diagnostics for the entire population is through probability-probability plots or P-P plots of the empirical distribution of the (PIT or HPD) statistic versus its distribution under the hypothesis that ; henceforth, we will refer to the latter Uniform(0,1) distribution as the “theoretical” distribution of PIT or HPD. An ideal P-P plot has all points close to the identity line where the “empirical” and “theoretical” distributions are the same. Note that HPD P-P plots, in particular, are valuable calibration tools if our goal is to calibrate the estimated densities so that the computed predictive regions have the right coverage.

Figure 3 (c)-(d) contains example P-P plots for the PIT and HPD values of photo- estimators, where the empirical (observed) distribution is close to the theoretical (ideal) distribution.

Figure 2: Schematic diagram of the construction of the Probability Integral Transform (PIT, left) and the Highest Probability Density (HPD, right) values for the estimated density at , where is the response at . In the plot to the right, the highlighted segments on the -axis form the so-called highest density region (HDR) of . The PIT and HPD values correspond to the area of the tail versus highest density region, respectively, of the estimate; here indicated by the blue versus green shaded areas.

4 Applications in Astronomy

We demonstrate8 the breadth of our CDE methods in three different astronomical use cases.

  1. Photo- estimation: univariate response with multivariate input. This is the standard prediction setting in which all four methods apply. In Section 4.1, we apply the methods to the Teddy photometric redshift data (Beck et al., 2017), and illustrate the need for loss functions to properly assess PDF estimates of redshift given photometric colors .

  2. Likelihood-free cosmological inference: multivariate response. For multiple response components we want to model the often complicated dependencies between these components; this is in contrast to approaches which model each component separately, implicitly introducing an assumption of conditional independence. In Section 4.2, we use an example of LFI for simulated weak lensing shear data to show how NNKCDE and RFCDE can capture more challenging bivariate distributions with curved structures; in this toy example represents cosmological parameters in the CDM-model, and represents (coarsely binned) weak lensing shear correlation functions.

  3. Spec- estimation: functional input. Standard prediction methods, such as random forests, do not typically fare well with functional features, simply treating them as unordered vectorial data and ignoring the functional structure. However, often there are substantial benefits to explicitly taking advantage of that structure as in fRFCDE. In Section 4.3, we compare the performance of a vectorial implementation of RFCDE with fRFCDE and FlexCode-Series for a spectroscopic sample from the Sloan Digital Sky Survey (SDSS; Alam et al. 2015). The input is here a high-resolution spectrum of a galaxy, and the response is the galaxy’s redshift .

Throughout this section we will report CDE loss mean and standard error for each method, using Gaussian confidence intervals for performance comparison. Since the CDE loss is an empirical mean and the test size is reasonably large in all the examples, the validity of this approximation is guaranteed by the central limit theorem.

4.1 Photo- estimation: Univariate response with multivariate input

Here we estimate photo- posterior PDFs using representative training and test data (Samples A and B, respectively) from the Teddy catalog by Beck et al. (2017).9 These data include 74309 training and 74557 test observations. Each observation has five features: the magnitude of the -band and the pairwise differences or “colors” , , , and .

Among the chief sources of uncertainty affecting photo-s estimated by ML techniques are the incompleteness and nonrepresentativity of training sets, defined by the mismatch in the distributions of training and test data in and , which may be extreme to the point of not guaranteeing mutual coverage. Realistically modeling incompleteness is highly challenging, requiring both simulations of SEDs and of the observational conditions of a given survey, which is outside the scope of this work. Accounting for redshift incompletenss is not a lost cause and may be accomplished by extrapolating outside of the training range by abandoning standard instance-based ML algorithms (see e.g., Leistedt & Hogg (2016)). Certain types of selection bias, known as covariate shift, can also be corrected by importance weights in the CDE loss (Izbicki et al., 2017; Freeman et al., 2017); see Appendix D for details and code. However, for simplicity, in this paper we consider only representative training sets with no disparity in color-space coverage, putting this demonstration on equal footing with all previous comparisons of photo- PDF methods.

We fit NNKCDE, RFCDE, DeepCDE, and FlexZBoost to the data. In addition we fit an RFCDE model, “RFCDE-Limited,” restricted to the first three of the five features, as a toy model which fails to extract some information from the features, allowing us to showcase the difference between PIT or HPD diagnostics and the CDE loss function. For DeepCDE we use a three-layer deep neural network with linear layers and reLu activations with neurons per layers, trained for epochs with Adam (Kingma & Ba, 2014). As our goal is to showcase its applicability we do not optimize the neural architecture (e.g., number of layers, number of neurons per layer, activation functions) nor the learning parameters (i.e., number of epochs, learning rate, momentum). To illustrate our validation methods, we also include the marginal distribution as an estimate of individual photo- distributions . This estimate will be the same regardless of .

Figure 3: (a)-(b) Histograms of PIT and HPD values for Teddy photo- data. Both versions of RFCDE as well as the marginal distribution have a uniform distribution indicating that the conditional density estimates are well-calibrated on average. FlexZBoost exhibits an overrepresentation of 0 and 1 values which indicates overly narrow CDEs. (c)-(d) Probability-Probability (P-P) plots of HPD and PIT values. Both sets of values are computed over data in the test set and their empirical distribution is plotted against the uniform distribution; i.e., the “theoretical distribution” of the HPD and PIT values when . If the estimated CDEs are well calibrated, the empirical and theoretical distributions should coincide and all points should be close to the identity line. As in the top panel, the P-P plots indicate a good fit for all methods including the clearly misspecified “Marginal” model.
Method CDE Loss SE Storage (single CDE) Marginal -3.192 0.007 200 floats (1.6 KB) RFCDE-Limited -11.038 0.047 200 floats (1.6 KB) NNKCDE -12.139 0.043 200 floats (1.6 KB) FlexZBoost -12.272 0.044 30 coeffs (0.24 KB) DeepCDE -12.821 0.04 31 coeffs (0.25 KB) RFCDE -13.108 0.052 200 floats (1.6 KB)
Figure 4: Method comparison via the CDE loss of Equation 5 with estimated standard error (SE) and the storage space for each galaxy’s photo- CDE. As FlexZBoost and DeepCDE are basis expansion methods, we need only to store the estimated coefficients for a lossless compression of the CDEs; the other CDEs are discretized to 200 bins.

Figure 4 presents the estimated CDE loss of Equation 5 with estimated standard error (SE), as well as the storage space for each galaxy’s photo- CDE. Figure 3 shows that the different models, including the clearly misspecified “Marginal” model, achieve comparable performance on goodness-of-fit diagnostics. However, Figure 4 shows the discriminatory power of the CDE loss function, which distinguishes the methods from one another. This emphasizes the need for method comparison through loss functions in addition to goodness-of-fit diagnostics. We note that the ranking in CDE loss also correlates roughly with the quality of the point estimates, as shown in Appendix B. The LSST-DESC PZ DC1 paper (Schmidt et al., 2019) draws similar conclusions from a comprehensive photo- code comparison where the “Marginal” model (there referred to as trainZ photo- PDF estimator, the experimental control) outperformed all codes when using traditional metrics for assessing photo- PDF acccuracy. Indeed, of the metrics considered in DC1, the CDE loss was the only metric that could appropriately penalize the pathological trainZ.

4.2 Likelihood-free cosmological inference: Multivariate response

To showcase the ability to target joint distributions, we apply NNKCDE and RFCDE to the problem of estimating multivariate posteriors of the cosmological parameters in a likelihood-free setting via ABC-CDE (Izbicki et al., 2019).

ABC is an approach to parameter inference in settings where the likelihood is not tractable but we have a forward model that can simulate data under fixed parameter settings . The simplest form of ABC is the ABC rejection sampling algorithm, where a set of parameters is first drawn from a prior distribution. The simulated data is accepted with tolerance if for some distance metric (e.g., the Euclidean distance) that measures the discrepancy between the simulated data and observed data . The outcome of the ABC rejection algorithm for small enough is a sample of parameter values approximately distributed according to the desired posterior distribution .

The basic idea of ABC-CDE is to improve the ABC estimate — and hence reduce the number of required simulations — by using the ABC sample/output as input to a CDE method tuned with the CDE loss restricted to a neighborhood of defined by the tolerance (Izbicki et al. (2019), Equation 3). Hence, in ABC-CDE, our CDE method can be seen as a post-adjustment method: it returns an estimate which we evaluate at the point to obtain a more accurate approximation of . This could also be beneficial in an active learning setting where the posterior distribution is used to identify relevant regions of the parameter space (e.g., (Lueckmann et al., 2019; Papamakarios et al., 2019; Alsing et al., 2019)).

In this example, we consider the problem of cosmological parameter inference via cosmic shear, caused by weak gravitational lensing inducing a distortion in images of distant galaxies. The size and direction of the distortion is directly related to the size and shape of the matter distribution along the line of sight, which varies across the universe. We use shear correlation functions to constrain the dark matter density and matter power spectrum normalization parameters of the CDM cosmological model, which predicts the properties and evolution of the large scale structure of the matter distribution. For further background see Hoekstra & Jain (2008), Munshi et al. (2008) and Mandelbaum (2017). Here we use the GalSim10 toolkit (Rowe et al., 2015) to generate simplified galaxy shears distributed according to a Gaussian random field determined by . The binned shear correlation functions serve as our input data or summary statistics . For the inference, we assume uniform priors and and fix , , .

Figure 5: Cosmological parameter inference in an LFI setting with simulated weak lensing data. The top row shows the estimated bivariate distribution of and for ABC rejection sampling at different acceptance rates (20%, 50%, and 100%). The middle row shows the estimated posterior densities after applying NNKCDE to the ABC sample, and the bottom row when applying RFCDE. Both NNKCDE and RFCDE tuned with a CDE loss improve on ABC across acceptance rate levels, i.e., they provide approximate posteriors that are more concentrated around the true observed parameter values and the degeneracy curve on which the data are indistinguishable (shown here as a black diamond and dashed line, respectively). We even see some structure at an ABC acceptance rate of 1 (right column); that is, at an ABC threshold of for which the entire sample is accepted by ABC and passed to our CDE code. HPD values were not included as they fail to distinguish between conditional and marginal distribution, as mentioned above.

The top row of Figure 5 shows the estimated bivariate posterior distribution of from ABC rejection sampling only, at varying acceptance rates (20%, 50%, and 100%). An acceptance rate of 100% just returns the (uniform) ABC prior distribution, whereas the ABC posteriors for smaller acceptance rates (that is, smaller values of ) concentrate around the parameter degeneracy curve (shown as a dashed line) on which the data are indistinguishable. The second and third rows show the estimated posteriors when we, respectively, improve the initial ABC estimate by applying NNKCDE and RFCDE tuned with our CDE surrogate loss. A result that is apparent from the figure is that NNKCDE and RFCDE fitted with a CDE loss are able to capture the degeneracy curve at a larger acceptance rate, that is, for a smaller number of simulations, than when using ABC only.

CDE Loss SE ()
Method \Acceptance Rate 20% 50% 100%
ABC -0.686 0.009 -0.392 0.004 -0.227 0.001
NNKCDE -1.652 0.022 -1.199 0.016 -0.844 0.010
RFCDE -4.129 0.063 -3.698 0.064 -2.817 0.055
Table 2: Performance of ABC, NNKCDE, RFCDE in LFI settings with simulated weak lensing data in terms of the surrogate CDE loss.

We also have a direct measure of performance via the CDE loss and can adapt to different types of data by leveraging different CDE codes. Table 2 shows how the methods compare in terms of the surrogate loss. While decreasing the acceptance rate benefits all methods, it is clear that CDE-based approaches have better performance in all cases.

4.3 Spec- estimation: Functional input

In this example, we compare CDE methods in the context of spectroscopic redshift prediction using 2812 high-resolution spectra from the Sloan Digital Sky Survey (SDSS) Data Release 6 (preprocessed with the cuts of Richards et al. 2009), corresponding to features of flux measurements at wavelengths. The high-resolution spectra can be seen as functional inputs on a continuum of wavelength values. Spectroscopic redshifts (or redshifts predicted from spectra ) tend to be both accurate and precise, so the density is well-approximated by a delta function at the true redshift. For the purposes of illustrating the use of the CDE codes, we define a noisified redshift , where are independent and identically distributed variables drawn from a normal distribution and is the true redshift of galaxy provided by SDSS. Thus the conditional density of this example is a Gaussian distribution with mean and variance 0.02.

We compare fRFCDE, a “Functional” adaptation of RFCDE, with a standard “Vector” implementation of RFCDE, which treats functional data as a vector. For completeness, we also compare against standard regression random forest combined with KDE, as well as FlexCode-Spec, an extension of FlexCode with a Spectral Series regression (Richards et al., 2009; Freeman et al., 2009; Lee & Izbicki, 2016) for determining the expansion coefficients. We train on 2000 galaxies and test on the remaining galaxies. The RFCDE and the fRFCDE trees are both trained with , , and bandwidths chosen by plug-in estimators. We use for the fRFCDE rate parameter of the Poisson process that defines the variable groupings.

Method Train Time (in sec) CDE Loss SE
Functional RFCDE 24.89 -3.38 0.155
Vector RFCDE 41.60 -2.52 0.104
Regression RF + KDE 50.63 -3.42 0.120
FlexCode-Spec 4017.93 -3.53 0.136
Table 3: Performance of RFCDE for functional data. We achieve both a lower CDE loss and computational time for Functional RFCDE(as compared to Vector RFCDE) by leveraging the functional nature of the data.

Table 3 contains the CDE loss and train time for both the vector-based and functional-based RFCDE models on the SDSS data, as well as for FlexCode-Spec. We obtain substantial gains when incorporating functional features both in terms of CDE loss and computational time. The computational gains are attributed to requiring fewer searches for each split point as the default value of is reduced. As anticipated, vector-based RFCDE underperforms with functional data due to the tree splits on CDE loss struggling to pick up signals in the data and returning almost the same conditional density regardless of the input. In contrast, splitting on mean squared error is an easier task and we include the results of regular regression random forest and KDE, which are comparable to the results of fRFCDE. Whereas random forests rely on variable selection (either individual variables as in vectorial RFCDE or grouped together as in fRFCDE) for dimensional reduction, FlexCode-Spec is based on the Spectral Series mechanism for dimension reduction that finds (potentially) nonlinear, sparse structure in the data distribution (Izbicki & Lee, 2016). Both types of dimension reduction are reasonable for spec- estimation, as particular wavelength locations or regions in the galaxy SED could carry information of the galaxy’s true redshift; RFCDE and fRFCDE are effective in finding such locations by variable selection. Spectral Series on the other hand are able to recover low-dimensional manifold structure in the entire data ensemble; as Richards et al. 2009 show, the main direction of variation in SDSS spectra is directly related to the spectroscopic redshift.

5 Conclusions

This paper presents statistical tools and software for uncertainty quantification in complex regression and parameter inference tasks. Given a set of features with associated response , our methods extend the usual point prediction task of classification and regression to nonparametric estimation of the entire (conditional) probability density . The described CDE methods are meant to handle a range of different data settings as outlined in Table 1. This paper includes examples of code usage in the contexts of photo- estimation, likelihood-free inference for cosmology analysis, and spec-z estimation. In addition, it provides tools for CDE method assessment and for choosing tuning parameters of the CDE methods in a principled way.

Our software includes four packages for CDE — NNKCDE, RFCDE, FlexCode and DeepCDE— each using a different machine learning algorithm, as well as a package for model assessment, cdetools. All packages are implemented in both Python and R, with no dependence on proprietary software, making them compatible with any operating system supporting either language (e.g., Windows, MacOS, Ubuntu, and any other Unix-based OS). The code is provided in publicly available Github repositories, documented using Python and R function documentation standards, and equipped with Travis CI11 automatic bug-tracking. Finally, our software makes uncertainty quantification straightforward for those used to standard open-source machine learning Python packages: NNKCDE, RFCDE, FlexCode share the sklearn API (with fit and predict methods), which makes our methods compatible with sklearn wrapper functions (for e.g., cross-validation and model ensembling). DeepCDE has implementations for both Tensorflow and Pytorch, two of the most widely used deep learning frameworks, and can easily be combined with most existing network architectures.


We would like to thank the two anonymous reviewers for their insightful comments that helped improve the manuscript. ND, TP and ABL were partially supported by the National Science Foundation under Grant Nos. DMS1521786. RI was supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (grant number 306943/2017-4) and Fundação de Amparo à Pesquisa do Estado de São Paulo (grants number 2017/03363-8 and 2019/11321-9). AIM acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. During the completion of this work, AIM was advised by David W. Hogg and supported in part by National Science Foundation grant AST-1517237.

Appendix A Examples of code usage

In this section we provide examples of code usage in Python and R. To save space, we have included the DeepCDE examples directly in the code release12 (Dalmasso & Pospisil, 2019) rather than in this section.

a.1 Nnkcde

{},codes=] import numpy as np import nnkcde

# Fit the model model = nnkcde.NNKCDE(), y_train)

# Tune parameters: bandwidth, number of neighbors k # Pick best loss on validation data k_choices = [5, 10, 100, 1000] bandwith_choices = [0.01, 0.05, 0.1] model.tune(x_validation, y_validation, k_choices, bandwidth_choices)

# Predict new densities on grid y_grid = np.linspace(y_min, y_max, n_grid) cde_test = model.predict(x_test, y_grid)  

{},codes=] library(NNKCDE)

# Fit the model model <- NNKCDE::NNKCDE$new(x_train, y_train)

# Tune parameters: bandwidth, number of neighbors k # Pick best loss on validation data k_choices <- c(5, 10, 100, 1000) bandwith_choices <- c(0.01, 0.05, 0.1) model$tune(x_validation, y_validation, k_grid = k_choices, h_grid = bandwidth_choices)

# Predict new densities on grid y_grid <- seq(y_min, y_max, length.out = n_grid) cde_test <- model$predict(x_test, y_grid)  

a.2 Rfcde

{},codes=] import numpy as np import rfcde

# Parameters n_trees = 1000 # Number of trees in the forest mtry = 4 # Number of variables to potentially split at in each node node_size = 20 # Smallest node size n_basis = 15 # Number of basis functions bandwidth = 0.2 # Kernel bandwith - used for prediction only

# Fit the model forest = rfcde.RFCDE(n_trees=n_trees, mtry=mtry, node_size=node_size, n_basis=n_basis) forest.train(x_train, y_train)

# Predict new densities on grid y_grid = np.linspace(0, 1, n_grid) cde_test = forest.predict(x_test, y_grid, bandwidth)

# Predict conditional means for CDE-optimized forest cond_mean_test = forest.predict_mean(x_test)

# Predict quantiles (i.e., quantile regression) for CDE-optimized forest alpha_quantile = .90 quant_test = forest.predict_quantile(x_test, alpha_quantile)  

{},codes=] library(RFCDE)

# Parameters n_trees <- 1000 # Number of trees in the forest mtry <- 4 # Number of variables to potentially split at in each node node_size <- 20 # Smallest node size n_basis <- 15 # Number of basis functions bandwidth <- 0.2 # Kernel bandwith - used for prediction only

# Fit the model forest <- RFCDE::RFCDE(x_train, y_train, n_trees = n_trees, mtry = mtry, node_size = node_size, n_basis = n_basis)

# Predict new densities on grid y_grid <- seq(0, 1, length.out = n_grid) cde_test <- predict(forest, x_test, y_grid, response = ’CDE’, bandwidth = bandwidth)

# Predict conditional means for CDE-optimized forest cond_mean_test <- predict(forest, x_test, y_grid, response = ’mean’, bandwidth = bandwidth)

# Predict quantiles (i.e., quantile regression) for CDE-optimized forest alpha_quantile <- .90 quant_test <- predict(forest, x_test, y_grid, response = ’quantile’, quantile = alpha_quantile, bandwidth = bandwidth)  


{},codes=] import numpy as np import rfcde

# Parameters n_trees = 1000 # Number of trees in the forest mtry = 4 # Number of variables to potentially split at in each node node_size = 20 # Smallest node size n_basis = 15 # Number of basis functions bandwidth = 0.2 # Kernel bandwith - used for prediction only lambda_param = 10 # Poisson Process parameter

# Fit the model functional_forest = rfcde.RFCDE(n_trees=n_trees, mtry=mtry, node_size=node_size, n_basis=n_basis) functional_forest.train(x_train, y_train, flamba=lambda_param)

# … Same as RFCDE for prediction …  

{},codes=] library(RFCDE)

# Parameters n_trees <- 1000 # Number of trees in the forest mtry <- 4 # Number of variables to potentially split at in each node node_size <- 20 # Smallest node size n_basis <- 15 # Number of basis functions bandwidth <- 0.2 # Kernel bandwith - used for prediction only lambda_param <- 10 # Poisson Process parameter

# Fit the model functional_forest <- RFCDE::RFCDE(x_train, y_train, n_trees = n_trees, mtry = mtry, node_size = node_size, n_basis = n_basis, flambda = lambda_param)

# … Same as RFCDE for prediction …  

a.3 FlexCode

{},codes=] import numpy as np import flexcode

# Select regression method from flexcode.regression_models import NN

# Parameters basis_system = "cosine" # Basis system max_basis = 31 # Maximum number of basis. If the model is not tuned, # max_basis is set as number of basis

# Regression Parameters # If a list is passed for any parameter automatic 5-fold CV is used to # determine the best parameter combination. params = {"k": [5, 10, 15, 20]} # A dictionary with method-specific regression parameters.

# Parameterize model model = flexcode.FlexCodeModel(NN, max_basis, basis_system, regression_params=params)

# Fit model - this will also choose the optimal number of neighbors ‘k‘, y_train)

# Tune model - Select the best number of basis model.tune(x_validation, y_validation)

# Predict new densities on grid cde_test, y_grid = model.predict(x_test, n_grid=n_grid)  

{},codes=] library(FlexCoDE)

# Parameters basis_system <- "Cosine" # Basis system max_basis <- 31 # Maximum number of basis.

# Fit and tune FlexCode via KNN regression using 4 cores fit <- FlexCoDE::fitFlexCoDE(x_train, y_train, x_validation, y_validation, nIMax = max_basis, regressionFunction.extra = list(nCores = 4), regressionFunction = FlexCoDE::regressionFunction.NN, system = basis_system)

# Predict new densities on a grid n_points_grid <- 500 # Number of points in the grid cde_test <- predict(fit, x_test, B = n_points_grid)

## Shortcut for FlexZBoost

# Fit and tune FlexZBoost fit <- FlexCoDE::FlexZBoost(x_train, y_train, x_validation, y_validation, nIMax = max_basis, system = basis_system)

# Predict new densities on a grid n_points_grid <- 500 # Number of points in the grid cde_test <- predict(fit, x_test, B = n_points_grid)  

a.4 cdetools

{},codes=] import cdetools

cde_test # numpy matrix of conditional density evaluations on a grid # - each row is an observation, each column is a grid point y_grid # the grid at which cde_test is evaluated y_true # the observed y values

# Calculate the cde_loss cde_loss(cde_test, y_grid, y_test)

# Calculate PIT values cdf_coverage(cde_test, y_grid, y_test)

# Calculate HPD values hpd_coverage(cde_test, y_grid, y_test)  

{},codes=] library(cde_tools)

cde_test # matrix of conditional density evaluations on a grid # - each row is an observation, each column is a grid point y_grid # the grid at which cde_test is evaluated y_true # the observed y values

# Calculate the cde_loss cde_loss(cde_test, y_grid, y_test)

# Calculate PIT values cdf_coverage(cde_test, y_grid, y_test)

# Calculate HPD values hpd_coverage(cde_test, y_grid, y_test)  

Appendix B Point Estimate Photo-Z Redshifts (Univariate Response with Multivariate Data)

Plots of point estimates can sometimes be a useful qualitative diagnostic of photo- code performance in that they may quickly identify a bad model. For example, the marginal model would always lead to the same point estimate and hence a horizontal line in the plot. Figure 6 shows point photo- predictions versus the true redshift for the photo- codes in Example 4.1.

Figure 6: Predicted photometric versus true spectroscopic redshift. Point estimate photo-s are derived from the CDEs by taking the mode of each density. Predictions from the marginal, not shown here, always predicts the mode of the marginal distribution, which here is .

The second moment of the distribution and outlier fraction (OLF) are also commonly used diagnostic tools for photo- estimators. Table 4 lists the , and OLF point estimate metrics from Dahlen et al. (2013). These metrics are defined as follows:


where indicates the spectroscopic redshift and is the difference between the predicted and spectroscopic redshift. According to the point estimate metrics, RFCDE and DeepCDE perform the best on the Teddy data, followed by FlexZBoost and NNKCDE. These results are consistent with the CDE loss rankings in Figure 4.

Method OLF (%)
Marginal 0.0049 0.0740 1.018
RFCDE-Limited 0.0010 0.0205 0.604
NNKCDE 0.0008 0.0180 0.335
FlexZBoost 0.0008 0.0169 0.362
DeepCDE 0.0007 0.0164 0.359
RFCDE 0.0007 0.0169 0.345
Table 4: Performance in photometric redshift predictions for different methods.

Appendix C Loss Function and Model Assumption Equivalence

In this section we connect the DeepCDE loss in Equation 4 to the FlexCode model setup in Equation 2.

Lemma 1.

Let be the coefficients of the FlexCode basis expansion in Equation 2. Minimizing the CDE loss in Equation 4 is equivalent to minimizing the mean squared errors of the basis expansion coefficients, i.e.,


We prove the statement by showing that the two minimization problems are equivalent.

First, considering the LHS of Equation C1, we have that:

where the second equality follows from the fact that the set is part of an orthonormal basis, that is,

Next, the RHS of Equation C1 reduces to:

In DeepCDE, we implement an especially simple expression for the loss function by inserting the orthogonal basis expansion (Equation 3) into the CDE loss (Equation 4) to yield


The last expression represents our empirical CDE loss on validation data and is easy to compute.

Appendix D Accounting for Covariate Shift in Photometric Redshift Estimation via Importance Weights

In the photo- application of Section 4.1 we have for simplicity assumed the distributions of the spectroscopic and photometric data to be the same. In this section we present a solution by Izbicki et al. (2017) and Freeman et al. (2017) for a setting where the spectroscopically labeled galaxies and the photometric “target” galaxies of interest lie in different regions of the color space. Under so-called covariate shift the the (marginal) distributions and of the spectroscopic and photometric data, respectively, are different — but their conditional distributions and are assumed to be the same. In other words:

Covariate shift arises from missing at random (MAR) bias, which for example occurs when the probability that a galaxy is labeled with a spectroscopic redshift depends only on its photometry and possibly other characteristics but not on its true redshift . In such a setting the goal is to evaluate the performance on the photometric “target” data of interest. Hence, in terms of loss functions, our goal is to minimize a CDE loss computed on the photometric, and not the spectroscopic, data as in


where the weight function will enforce a good fit in the regions of color or feature space where most of the photometric data reside. The challenge is how to compute the above loss for the photometric data of interest, as the true redshift , by construction, is only known for the spectroscopic set. A solution proposed in the above-mentioned papers, is to rewrite the loss by introducing so-called importance weights . An example of how to estimate such weights can be found in Izbicki et al. (2017), Sugiyama et al. (2007), Gretton et al. (2009) and Kremer et al. (2015). Once the importance weights have been estimated, one can rewrite an estimate of the loss in D1 as an expression that can be computed on held-out data without knowing the redshift of the (photometric) target data:


where and represent the spectroscopic and photometric data, respectively, in the held-out set. Furthermore, as noted in Izbicki et al. (2017) (Section 5), one can often drastically improve CDE performance by just restricting the spectroscopic training data to regions where the unlabeled photometric data reside. A simple procedure for identifying such training examples is to compute importance weights for all labeled examples, and then replace the examples in the training set with with new labeled examples that have . Effectively, this will create a spectroscopic training set that resides in the same part of the feature space as the unlabeled photometric test data (when that is possible).

Figure 7 (adapted from Izbicki et al. (2017)) shows the performance of three versions of the FlexCode-series. The difference between these approaches is in how they are tuned. The first result (“no reweighting”) corresponds to the standard version of FlexCode that does not take into account covariate shift and hence is tuned on held-out validation data according to Equation (5). The second estimator (“removing zeros”) is also tuned according to Equation (5), but on a held-out sample that has no training examples in regions with . For the last estimator (“reweighting”), tuning is performed according to Equation (D2). The plot shows that the simple approach of removing galaxies with zero estimated weights already improves the estimators, but reweighting the loss by can yield even better performance.

Figure 7: Performance of three versions of FlexCode-Series to the photometric data used in (Izbicki et al., 2017). Left: estimator that ignores selection bias; middle: removing and replacing galaxies with in the training set, but otherwise not including the importance weights in the computations; right: reweighting galaxies by . (Adapted from (Izbicki et al., 2017)).


  1. All of these methods, except for DeepCDE, have occurred in previously published or archived papers. Hence, we only briefly review the methods in this paper and instead focus on software usage, algorithmic aspects, and the settings under which each method can be applied.
  4. Available at
  5. (Python) and (R)
  8. Code for these examples is publicly available at
  9. Data available at


  1. Abadi M., et al., 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems,
  2. Abbott T. M. C., et al., 2019, Phys. Rev. Lett., 122, 171301
  3. Aghanim N., et al., 2018, arXiv preprint arXiv:1807.06209
  4. Alam S., et al., 2015, The Astrophysical Journal Supplement Series, 219, 12
  5. Alsing J., Wandelt B., Feeney S., 2018, Monthly Notices of the Royal Astronomical Society, 477, 2874
  6. Alsing J., Charnock T., Feeney S., Wandelt B., 2019, arXiv preprint arXiv:1903.00007
  7. Ball N. M., Brunner R. J., 2010, International Journal of Modern Physics D, 19, 1049
  8. Beaumont M. A., Zhang W., Balding D. J., 2002, Genetics, 162, 2025
  9. Beck R., et al., 2017, Monthly Notices of the Royal Astronomical Society, 468, 4323
  10. Bishop C. M., 1994, Technical report, Mixture density networks. Citeseer
  11. Breiman L., 2001, Machine learning, 45, 5
  12. Carrasco Kind M., Brunner R. J., 2013, Monthly Notices of the Royal Astronomical Society, 432, 1483
  13. Carrasco Kind M., Brunner R. J., 2014, Mon Not R Astron Soc, 441, 3550
  14. Chen T., Guestrin C., 2016, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp 785–794, doi:10.1145/2939672.2939785
  15. Chen Y., Gutmann M. U., 2019, in Chaudhuri K., Sugiyama M., eds, Proceedings of Machine Learning Research Vol. 89, Proceedings of Machine Learning Research. PMLR, pp 1584–1592,
  16. Cunha C. E., Lima M., Oyaizu H., Frieman J., Lin H., 2009, Monthly Notices of the Royal Astronomical Society, 396, 2379
  17. D’Isanto A., Polsterer K., 2018, Astronomy & Astrophysics, 609, A111
  18. Dahlen T., et al., 2013, The Astrophysical Journal, 775, 93
  19. Dalmasso N., Pospisil T., 2019, tpospisi/DeepCDE 0.1, doi:10.5281/zenodo.3364862,
  20. Dalmasso N., Pospisil T., Izbicki R., 2019, tpospisi/cdetools 0.0.1, doi:10.5281/zenodo.3364810,
  21. Epanechnikov V. A., 1969, Theory of Probability & Its Applications, 14, 153
  22. Fischer P., Dosovitskiy A., Brox T., 2015, in GCPR. , doi:10.1007/978-3-319-24947-6_30
  23. Frank E., Hall M., 2001, in De Raedt L., Flach P., eds, Machine Learning: ECML 2001. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 145–156
  24. Freeman P. E., Newman J. A., Lee A. B., Richards J. W., Schafer C. M., 2009, Monthly Notices of the Royal Astronomical Society, 398, 2012
  25. Freeman P. E., Izbicki R., Lee A. B., 2017, Monthly Notices of the Royal Astronomical Society, 468, 4556
  26. Genest C., Rivest L.-P., 2001, Statistics & Probability Letters, 53, 391
  27. Goodfellow I., Bengio Y., Courville A., 2016, Deep learning. MIT press
  28. Graff P., Feroz F., Hobson M. P., Lasenby A., 2014, Monthly Notices of the Royal Astronomical Society, 441, 1741
  29. Greenberg D. S., Nonnenmacher M., Macke J. H., 2019, in Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA. pp 2404–2414,
  30. Gretton A., et al., 2009, Dataset Shift in Machine Learning, 131-160 (2009)
  31. Harrison D., Sutton D., Carvalho P., Hobson M., 2015, Monthly Notices of the Royal Astronomical Society, 451, 2610
  32. Hoekstra H., Jain B., 2008, Annual Review of Nuclear and Particle Science, 58, 99
  33. Hyndman R. J., 1996, The American Statistician, 50, 120
  34. Hyndman R., Einbeck J., Wand M., 2018, hdrcde: Highest Density Regions and Conditional Density Estimation.
  35. Ivezić Ž., Connolly A. J., VanderPlas J. T., Gray A., 2014, Statistics, data mining, and machine learning in astronomy: a practical Python guide for the analysis of survey data. Vol. 1, Princeton University Press
  36. Ivezić Ž., et al., 2019, The Astrophysics Journal, 873, 111
  37. Izbicki R., Lee A. B., 2016, Journal of Computational and Graphical Statistics, 25, 1297
  38. Izbicki R., Lee A. B., 2017, Electronic Journal of Statistics, 11, 2800
  39. Izbicki R., Pospisil T., 2019, rizbicki/FlexCode v5.9-beta.3, doi:10.5281/zenodo.3366065,
  40. Izbicki R., Lee A., Schafer C., 2014, in Artificial Intelligence and Statistics. pp 420–429
  41. Izbicki R., Lee A. B., Freeman P. E., 2017, The Annals of Applied Statistics, 11, 698
  42. Izbicki R., Lee A. B., Pospisil T., 2019, Journal of Computational and Graphical Statistics, pp 1–20
  43. Joudaki S., et al., 2017, Monthly Notices of the Royal Astronomical Society, 474, 4894
  44. Juric M., et al., 2017, Data Products Definition Document,
  45. Kingma D. P., Ba J., 2014, CoRR, abs/1412.6980
  46. Krause E., et al., 2017, arXiv preprint arXiv:1706.09359
  47. Kremer J., Gieseke F., Pedersen K., Igel C., 2015, Astronomy and Computing, 193
  48. Krizhevsky A., Sutskever I., Hinton G. E., 2012, in Proceedings of the 25th International Conference on Neural Information Processing Systems - Volume 1. NIPS’12. Curran Associates Inc., USA, pp 1097–1105,
  49. LeCun Y., Bengio Y., Hinton G., 2015, nature, 521, 436
  50. Lee A. B., Izbicki R., 2016, Electronic Journal of Statistics, 10, 423
  51. Leistedt B., Hogg D., 2016, The Astrophysical Journal, 838
  52. Li J., Nott D., Fan Y., Sisson S., 2017, Computational Statistics & Data Analysis, 106, 77
  53. Lueckmann J.-M., Gonçalves P. J., Bassetto G., Öcal K., Nonnenmacher M., Macke J. H., 2017, in Proceedings of the 31st International Conference on Neural Information Processing Systems. NIPS’17. Curran Associates Inc., USA, pp 1289–1299,
  54. Lueckmann J.-M., Bassetto G., Karaletsos T., Macke J. H., 2019, in Symposium on Advances in Approximate Bayesian Inference. pp 32–53
  55. Malz A. I., Marshall P. J., DeRose J., Graham M. L., Schmidt S. J., Wechsler R., Collaboration) L. D. E. S., 2018, AJ, 156, 35
  56. Mandelbaum R., 2017, arXiv preprint:1710.03235
  57. Mandelbaum R., et al., 2008, Monthly Notices of the Royal Astronomical Society, 386, 781/806
  58. Marin J.-M., Raynal L., Pudlo P., Ribatet M., Robert C., 2016, Bioinformatics (Oxford, England), 35
  59. Meinshausen N., 2006, Journal of Machine Learning Research, 7, 983
  60. Munshi D., Valageas P., Van Waerbeke L., Heavens A., 2008, Physics Reports, 462, 67
  61. Nadaraya E. A., 1964, Theory of Probability & Its Applications, 9, 141
  62. Ntampaka M., et al., 2019, arXiv preprint arXiv:1902.10159
  63. Papamakarios G., Murray I., 2016, in Lee D. D., Sugiyama M., Luxburg U. V., Guyon I., Garnett R., eds, Advances in Neural Information Processing Systems 29. Curran Associates, Inc., pp 1028–1036
  64. Papamakarios G., Sterratt D., Murray I., 2019, in The 22nd International Conference on Artificial Intelligence and Statistics. pp 837–848
  65. Pasquet J., Bertin E., Treyer M., Arnouts S., Fouchez D., 2019, Astronomy & Astrophysics, 621, A26
  66. Paszke A., et al., 2019, in Wallach H., Larochelle H., Beygelzimer A., d’ Alché-Buc F., Fox E., Garnett R., eds, , Advances in Neural Information Processing Systems 32. Curran Associates, Inc., pp 8024–8035
  67. Pedregosa F., et al., 2011, J. Mach. Learn. Res., 12, 2825
  68. Pereira C. A. d. B., Stern J., 1999, Entropy, 1, 99
  69. Polsterer K. L., D’Isanto A., Gieseke F., 2016, Uncertain Photometric Redshifts
  70. Pospisil T., Dalmasso N., 2019a, tpospisi/NNKCDE 0.3, doi:10.5281/zenodo.3364858,
  71. Pospisil T., Dalmasso N., 2019b, tpospisi/RFCDE 0.3.2, doi:10.5281/zenodo.3364856,
  72. Pospisil T., Lee A. B., 2018, arXiv preprint arXiv:1804.05753
  73. Pospisil T., Lee A. B., 2019, (f)RFCDE: Random Forests for Conditional Density Estimation and Functional Data (arXiv:1906.07177)
  74. Pospisil T., Dalmasso N., Inacio M., 2019, tpospisi/FlexCode 0.1.5, doi:10.5281/zenodo.3364860,
  75. Ravikumar P., Lafferty J., Liu H., Wasserman L., 2009, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71, 1009
  76. Richards J. W., Freeman P. E., Lee A. B., Schafer C. M., 2009, Astrophysical Journal, 691, 32
  77. Rowe B., et al., 2015, Astronomy and Computing, 10, 121
  78. Schmidt S., Malz A., et al. 2019, (In Preparation) Under internal review by LSST-DESC
  79. Sheldon E. S., Cunha C. E., Mandelbaum R., Brinkmann J., Weaver B. A., 2012, The Astrophysical Journal Supplement Series, 201, 32
  80. Sisson S. A., Fan Y., Beaumont M., 2018, Handbook of Approximate Bayesian Computation. Chapman and Hall/CRC, doi:10.1201/9781315117195
  81. Sohn K., Lee H., Yan X., 2015, in Cortes C., Lawrence N. D., Lee D. D., Sugiyama M., Garnett R., eds, , Advances in Neural Information Processing Systems 28. Curran Associates, Inc., pp 3483–3491,
  82. Sugiyama M., Nakajima S., Kashima H., Bünau P. v., Kawanabe M., 2007, in Proceedings of the 20th International Conference on Neural Information Processing Systems. NIPS’07. Curran Associates Inc., USA, pp 1433–1440,
  83. Tang Y., Salakhutdinov R. R., 2013, in Burges C. J. C., Bottou L., Welling M., Ghahramani Z., Weinberger K. Q., eds, , Advances in Neural Information Processing Systems 26. Curran Associates, Inc., pp 530–538,
  84. Tibshirani R., 1996, Journal of the Royal Statistical Society: Series B (Methodological), 58, 267
  85. Vanderplas J., Connolly A., Ivezić Ž., Gray A., 2012, in Conference on Intelligent Data Understanding (CIDU). pp 47 –54, doi:10.1109/CIDU.2012.6382200
  86. Viironen K., et al., 2018, Astronomy & Astrophysics, 614, A129
  87. Watson G. S., 1964, Sankhyā: The Indian Journal of Statistics, Series A, pp 359–372
  88. Way M. J., Scargle J. D., Ali K. M., Srivastava A. N., 2012, Advances in Machine Learning and Data Mining for Astronomy. Chapman and Hall/CRC
  89. Wittman D., 2009, The Astrophysical Journal Letters, 700, L174
  90. van Uitert E., et al., 2018, Monthly Notices of the Royal Astronomical Society, 476, 4662
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description