{'fraction_non_alphanumeric': 0.04047035636286156, 'fraction_numerical': 0.014195919671869174, 'mean_word_length': 4.7117973682208865, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
LinCDE: Conditional Density Estimation via Lindsey's Method

Conditional density estimation is a fundamental problem in statistics, with scientific and practical applications in biology, economics, finance and environmental studies, to name a few. In this paper, we propose a conditional density estimator based on gradient boosting and Lindsey's method (LinCDE). LinCDE admits flexible modeling of the density family and can capture distributional characteristics like modality and shape. In particular, when suitably parametrized, LinCDE will produce smooth and non-negative density estimates. Furthermore, like boosted regression trees, LinCDE does automatic feature selection. We demonstrate LinCDE's efficacy through extensive simulations and three real data examples.


Introduction
In statistics, a fundamental problem is characterizing how a response depends on a set of covariates. Numerous methods have been developed for estimating the mean response conditioning on the covariates-the so-called regression problem. However, the conditional mean may not always be sufficient in practice, and various distributional characteristics or even the full conditional distribution are called for, such as the mean-variance analysis of portfolios (Markowitz, 1959), the bimodality of gene expression distributions (DeSantis et al., 2014;Moody et al., 2019), and the peak patterns of galaxy redshift densities (Ball et al., 2008). Conditional distributions can be used for constructing prediction intervals, downstream analysis, visualization, and interpretation (Arnold et al., 1999). Therefore, it is worthwhile to take a step forward from the conditional mean to the conditional distribution.
There are several difficulties in estimating conditional distributions. First, distribution estimation is more complicated than mean estimation regardless of the conditioning. Second, as with conditional mean estimation, conditioning on a potentially large number of covariates suffers from the curse of dimensionality. When only a small subset of the covariates are relevant, proper variable selection is necessary to mitigate overfitting, reduce computational burden, and identify covariates that may be of interest to the practitioners.
In this paper, we develop a tree-boosted conditional density estimator based on Lindsey's method, which we call LinCDE (pronounced "linseed") boosting. LinCDE boosting is built on the base learner LinCDE tree. A LinCDE tree partitions the covariate space into subregions with homogeneous conditional distributions, estimates a local unconditional density in each subregion, and aggregates the unconditional densities to form the final conditional estimator. LinCDE boosting combines LinCDE trees to form a strong ensemble learner.
LinCDE boosting possesses several desirable properties. LinCDE boosting provides a flexible modeling paradigm and is capable of capturing distributional properties such as heteroscedasticity and multimodality. LinCDE boosting also inherits the advantages of tree and boosting methods, and in particular, LinCDE boosting is able to detect influential covariates. Furthermore, the conditional density estimates are automatically non-negative and smooth, and other useful statistics such as conditional quantiles or conditional cumulant distribution functions (CDFs) can be obtained in a straightforward way from the LinCDE boosting estimates.
The organization of the paper is as follows. In Section 2, we formulate the problem and discuss related work. We develop LinCDE in three steps: • In Section 3, we describe Lindsey's method for (marginal) density estimation.
• In Section 4, we introduce LinCDE trees for conditional density estimation, which combine Lindsey's method with recursive partitioning.
• In Section 5, we develop a boosted ensemble model using LinCDE trees.
In Section 6, we discuss two optional but helpful preprocessing steps-response transformation and conditional mean centering. In Section 7, we evaluate the performance of LinCDE boosting on simulated data sets. In Section 8, we apply LinCDE boosting to three real data sets. We conclude the paper with discussions in Section 9 and provide links to the R software and instructions in Section 10. All proofs are deferred to the appendix.

Background
In this section, we formulate the conditional density estimation problem and discuss related work.

Problem Formulation
Let y ∈ R be a continuous response. 1 Let x be a d-dimensional covariate vector and x (j) be its j-th coordinate. We assume the covariates are generated from an unknown underlying distribution f x (x), and the response given the covariates are sampled from an unknown conditional density f y|x (y | x). The model is summarized as ind.
∼ f y|x . (1) 1. The paper will focus on univariate responses, and the generalization to multivariate responses is straightforward and discussed in Section 9.
We observe n data pairs {(x i , y i )} and aim to estimate the conditional density f y|x (y | x).

Literature
In general, there are three ways to characterize a conditional distribution: conditional density, conditional quantile, and conditional CDF. We categorize related works on conditional density, quantile, and CDF estimation according to the methodology and provide a review below. We conclude by discussing two desired properties of conditional density estimation and how our work complements the existing literature on the two points. A line of study estimates the conditional distribution by localizing unconditional distribution estimators. Localization methods weight observations according to the distances between their covariates and those at the target point, and solve the unconditional distribution estimation problem based on the weighted sample. For conditional density, Fan et al. (1996) obtain conditional density estimates by local polynomial regression. For conditional quantile, Chaudhuri (1991a,b) partitions the covariate space into bins and fit a quantile model in each bin separately, and Yu and Jones (1998) tackle the conditional quantile estimation via local quantile loss minimization. For conditional CDF, Stone (1977) proposes a weighted sum of indicator functions, and Hall et al. (1999) consider a local logistic regression and a locally adjusted Nadaraya-Watson estimator. Localization methods enable systematic extensions of any unconditional estimator. Nevertheless, the weights usually treat covariates as equally important, and variable selection is typically not accommodated. This leaves the methods vulnerable to the curse of dimensionality.
As discussed by Stone (1991a,b), another approach making use of unconditional methods, first obtains the joint distribution estimatef y,x (y, x) and the covariate distribution estimatef x (x), and then followsf y|x (y | x) =f y,x (y, x)/f x (x) (2) to derive the conditional density. Nevertheless, the joint distribution estimation is also challenging, if not more. Arnold et al. (1999) point out except for special cases like multivariate Gaussian, the estimation of a bivariate joint distribution in a certain exponential form is onerous due to the normalizing constant. Moreover, the approach is inefficient both statistically and computationally if the conditional distribution is comparatively simpler than the joint and covariate distributions; as an example, when the response is independent of the covariates. A different thread directly models the dependence of the response's distribution on the covariates by a linear combination of a finite or infinite number of bases.
• For conditional density, Stone (1991, 1992), Stone (1991b), Mâacsse and Truong (1999), and Barron and Sheu (1991) study the conditional logspline density model: modeling log(f y|x (y | x)) by tensor products of splines or trigonometric series and maximizing the conditional log-likelihood to estimate the parameters. The method is also known as entropy maximization subject to empirical constraints. In addition, Sugiyama et al. (2010) model the conditional density in reproducing kernel Hilbert spaces and estimate the loadings by the unconstrained least-squares importance fitting, and Izbicki and Lee (2016) expand the conditional density in the eigen-functions of a kernel-based operator which adapts to the intrinsic dimension of the covariates. 2 • For conditional quantiles, Koenker and Bassett Jr (1978) formulate conditional quantiles as linear functions of covariates and minimize quantile losses to estimate parameters. Koenker et al. (1994) explore quantile smoothing splines minimizing for a single covariate, and He et al. (1998) extend the approach to the bi-variate setting, i.e., two covariates. Li et al. (2007) propose kernel quantile regression (KQR) considering quantile regression in reproducing Hilbert kernel spaces. Belloni et al. (2019) approximate the conditional quantile function by a growing number of bases as the sample size increases.
• For conditional CDF, Foresi and Peracchi (1995) and Chernozhukov et al. (2013) consider distribution regression: estimating a sequence of conditional logit models over a grid of values of the response variable. Belloni et al. (2019) further extend the method to the high-dimensional-sparse-model setting.
The performance of parametric models depends on the selected bases or kernels. Covariatespecific bases or kernels require a lot of tuning, and the bases or kernels treat covariates equally, making the approaches less powerful in the presence of many nuisance covariates. More recently, tree-based estimators arose in conditional distribution estimation. The overall idea is partitioning the covariate space recursively and fitting an unconditional model at each terminal node. For conditional quantiles, Chaudhuri and Loh (2002) investigate a tree-structured quantile regression. Nevertheless, the estimation of different quantiles requires separate quantile loss minimization, which complicates the full conditional distribution calculation. Meinshausen (2006) proposes the quantile regression forest (QRF) that computes all quantiles simultaneously. QRF first builds a standard random forest, then estimates the conditional CDF by a weighted distribution of the observed responses, and finally inverts the CDF to quantiles. For conditional CDF, Friedman (2019) proposes distribution boosting (DB). DB relies on Friedman's contrast trees-a method to detect the lack-of-fit regions of any conditional distribution estimator. DB estimates the conditional distribution by iteratively transforming the conditional distribution estimator and correcting the errors uncovered by contrast trees.
Beyond trees, neural network-based conditional distribution estimators have also been developed. For conditional quantile, standard neural networks with the quantile losses serve the goal. For conditional density, Bishop (1994Bishop ( , 2006 introduces the mixture density networks (MDNs) which model the conditional density as a mixture of Gaussian distributions. Neural network-based methods can theoretically approximate any conditional distributions well but are computationally heavy and lack interpretation.
In this paper, we propose LinCDE boosting which complements existing works by producing smooth density curves and performing feature selection.
• Compared to quantiles and CDF, smooth density curves are more suitable for the following reasons.
2. We remark that the approach in (Izbicki and Lee, 2016) adapts to the low intrinsic dimensionality of covariates but not the sparse dependency of the response on the covariates.
-Smooth density curves can give valuable indications of distributional characteristics such as skewness and multimodality. Visualization of smooth density curves is comprehensible to practitioners and can yield self-evident conclusions or point the way to further analysis. -Conditional densities can be used to compute class-posterior probabilities using the Bayes' rule. The class-posterior probabilities can be further used for classification. -Densities can be used to detect outliers: if an observation lies in a very lowdensity region, the data point is likely to be an outlier.
LinCDE boosting generates smooth density curves, while directly transforming estimates of conditional quantiles or CDF to conditional density estimates usually produces bumpy results.
• LinCDE trees and LinCDE boosting calculate feature importances and automatically focus on influential features in estimation. In contrast, methods like localization and modeling conditional densities in a linear space or RKHS often treat covariates equally, and covariate-specific neighborhoods, bases, and kernels may require a lot of tuning. Therefore, LinCDE boosting is more scalable to a large number of features and less sensitive to the presence of nuisance covariates.

Lindsey's Method
In this section, we first introduce the density estimation problem -an intermediate step towards the conditional density estimation. We then discuss how to solve the density estimation problem by Lindsey's method (Lindsey, 1974)-a stepping stone of LinCDE. Lindsey's method cleverly avoids the normalizing issue by discretization and solves the problem by fitting a simple Poisson regression. It can be thought of as a method for fitting a smooth histogram with a large number of bins. We consider the density family where κ : R → R is some carrying density, and g(y) is known as a tilting function. The idea is that κ(y) is known or assumed (such as Gaussian or uniform), and g(y) is represented by a model. We represent g as a linear expansion where z(y) is a basis of k smooth functions. As a simple example, if we use standard Gaussian as the carrying measure and choose z(y) = (y, y 2 ), the resulting density family corresponds to all possible Gaussian distributions. More generally we use a basis of natural cubic splines in z(y) with knots spread over the domain of y, to achieve a flexible representation (Wahba, 1975(Wahba, , 1990. Our goal is to find the density that maximizes the log-likelihood max β,β 0 1 n n i=1 log(κ(y i )) + z(y i ) β + β 0 , s.t. κ(y)e z(y) β+β 0 dy = 1.
The constrained optimization problem (5) can be simplified to the unconstrained counterpart below by the method of Lagrange multipliers (Silverman, 1986), where the multiplier turns out to be one. The optimization problem (6) is difficult since the integral κ(y)e z(y) β+β 0 dy is generally unavailable in closed form. One way to avoid the integral is by discretization-the key idea underlying Lindsey's method. One divides the response range into B equal bins of width ∆ with mid-points y b . The integral is approximated by the finite sum As for the first part of (6), one replaces y i by its bin midpoint y b(i) and groups the observations, where b(i) denotes the bin that the i-th response falls in, and n b represents the number of samples in bin b. Combining the above two parts, the Lagrangian function with response discretization takes the form The objective function (7) is equivalent to that of a Poisson regression with B observations {n b } 1≤b≤B and mean parameters µ b ∝ κ(y b )e z(y b ) β . Therefore, Lindsey's method estimates the coefficient β by fitting the Poisson regression with predictors z(y) and offset log (κ(y b )). The normalizing constant β 0 in (7) is absorbed in the Poisson regression's intercept. Despite the discretization error, Lindsey's estimates are consistent, asymptotically normal, and remarkably efficient (Moschopoulos and Staniswalis, 1994;Efron, 2004). We will demonstrate the efficacy of Lindsey's method in two examples at the end of this section. The number of bins B balances the statistical performance and the computational complexity of Lindsey's method: as B increases, the discretized objective (7) approaches the original target (6), and the resulting estimator converges to the true likelihood maximizer; on the other hand, the computations increase linearly in B. This can become a factor later when we fit many of these Poisson models repeatedly.
The relationship with a histogram becomes clear now, as well. We could use the counts in the B bins to form a density estimate, but this would be very jumpy. Typically we would control this by reducing the number of bins. Lindsey's method finesses this by having B large, but controlling the smoothness of the bin means via the k B basis functions and associated coefficients.
To control the model complexity further and avoid numeric instability, we add a regularization term to (6). For example, we can penalize deviations from Gaussian distributions via the regularizer (Wahba, 1977;Silverman, 1982Silverman, , 1986 The penalty measures the roughness of the tilting function and is zero if and only if the tilting function's exponent is a quadratic polynomial, i.e., a Gaussian distribution. We also attach a hyper-parameter λ to trade-off the objective (6) and the penalty (8), and tune λ to achieve the best performance on validation data sets. 4 It is convenient to tailor the spline basis functions to the penalty (8). Note that, for arbitrary bases z(y), the penalty (8) where Ω jl = z j (y)z l (y)dy. We transform our splines so that the associated Ω = diag(ω 1 , . . . , ω k ) is diagonal and the penalty reduces to a weighted ridge penalty k j=1 ω j β 2 j .
(details in Appendix B). Figure 1 depicts an example of the transformed spline bases (in increasing order of ω j ) and the corresponding smoothed versions. Among the transformed bases, the linear and quadratic components (the first and the second bases in Figure 1) are not shrunk by the roughness penalty (Claim 1), and higher-complexity splines are more heavily penalized (Hastie et al., 2009, chap. 5, for example).
Claim 1 Assume u ∈ R k and Ωu = 0. Then z(y) u is a linear or quadratic function of y.
We provide examples of Lindsey's performance in estimating bimodal and skewed distributions in Appendix G.

LinCDE Trees
In this section, we extend the density estimation problem to the conditional density estimation problem. We introduce LinCDE trees combining Lindsey's method and recursive partitioning: LinCDE trees partition the covariate space and estimate a local unconditional density via Lindsey's method in each subregion. To begin, we restate the target density family (3) in the language of exponential families. We call z(y) the sufficient statistics and g(y) the natural parameter. We replace the normalizing constant β 0 by the negative cumulant generating function ψ(β) defined as e ψ(β) = κ(y)e z(y) β dy.
As a result, the density normalizing constraint of κ(y)e z(y) β−ψ(β) is automatically satisfied.
In the conditional density estimation problem, we consider the target family generalized from (3) The dependence of the response on the covariates is encoded in the parameter function β(x).
• If there is only k = 1 basis function z(y) and the parameter function satisfies β(x) = x θ, then the conditional density family (11) reduces to a generalized linear model.
• If we consider a standard Gaussian prior, basis function z(y) = y, and tree-structured β(x), then the conditional density family (11) reduces to a regression tree model. In fact, the optimization problem based on (11) is equivalent to finding the tree β(x) that minimizes the sum of squares n i=1 (y i − β(x i )) 2 .
Similar to the density estimation problem, we aim to find the member in the family (11) that maximizes the conditional log-likelihood with the ridge penalty (10) where R 0 denotes the full covariate space. If β(x) is constant, the problem (12) simplifies to the unconditional problem (5). The conditional density estimation problem (12) is more complicated than the unconditional version due to the covariates x: 1. Given a covariate configuration x, there is often at most one observation whose covariates take the value x, and it is infeasible to estimate the multi-dimensional natural parameter β(x) based on a single observation; 2. There may be a multitude of covariates, and only a few are influential. Proper variable selection or shrinkage is necessary to avoid serious overfitting.
One way to finesse these difficulties is to use trees (Breiman et al., 1984). We divide the covariate space into subregions with approximately homogeneous conditional distributions, and in each subregion, we estimate a density independent of the covariate values. We name the method "LinCDE trees". In response to the first difficulty, by conditioning on a subregion instead of a specific covariate value, we have more samples for local density estimation. In response to the second difficulty, trees perform internal feature selection, and are thus resistant to the inclusion of many irrelevant covariates. Moreover, the advantages of tree-based methods are automatically inherited, such as being tolerant of all types of covariates, computationally efficient, and easy to interpret.
Before we delve into the details, we again draw a connection between LinCDE trees and a naive binning approach-fitting a multinomial model using trees. The naive approach discretizes the response into multiple bins and predicts conditional cell probabilities through recursive partitioning. The normalized conditional cell probabilities serve as an approximation of the conditional densities, and the more bins used, the higher resolution the approximation is. The naive approach is able to detect subregions with homogeneous multinomial distributions. However, the estimates are bumpy, especially with a large number of bins. To stabilize the method, restrictions enforcing smoothness are required, and LinCDE trees realize the goal by modeling the density exponent using splines.
We now explain how LinCDE trees work. In standard tree algorithms, there are two major steps: • Splitting: partitioning the covariate space into subregions; • Fitting: performing estimation in each subregion. The estimator is usually obtained by maximizing a specific objective function. For example, in a regression tree with 2 loss, the estimator is the sample average; in a classification tree with misclassification error, the estimator is the majority's label.
The fitting step is a direct application of Lindsey's method in Section 3. In a subregion R, we treat the natural parameter functions as a constant vector and solve the density estimation problem via Lindsey's method. We denote the objective function value in region R with parameter β by (R; β), and letβ R := arg max β (R; β). Now for the splitting step. Similar to standard regression and classification trees, we proceed with a greedy algorithm and select the candidate split that improves the objective the most. Mathematically, starting from a region R, we maximize the improvement statistic where R s,L and R s,R are the regions on the left and right of the candidate split, respectively. Direct computation of the difference (13) requires running Lindsey's method twice for each candidate split s to obtainβ R s,L ,β R s,R , and the total computation time is prohibitive. Instead, we approximate the difference (13) by a simple quadratic term in Proposition 2, which can be computed much faster.
Proposition 2 (Improvement approximation for LinCDE trees) Let n R ,z R be the sample size and average sufficient statistics in a region R. Assume that ∇ 2 ψ(β R ) + 2λΩ is invertible, then for a candidate split s, Proposition 2 writes the difference (13) as a quadratic form plus a higher-order residual term. If z(y) = y, the model amounts to a regression tree, and the residual term is zero. For general z(y), when the average sufficient statisticsz R s,L ,z R s,R are similar, the residual term is of smaller order than the quadratic form and can thus be dropped; whenz R s,L ,z R s,R are considerably different, the residual term is not guaranteed to be small theoretically. However, we empirically demonstrate in Appendix C that at such splits, the quadratic form is still sufficiently close to the true log-likelihood difference. Based on this empirical evidence, we use the quadratic approximation to determine the optimal splits. The quadratic approximation suggested by Proposition 2 is the product of the squared difference between the average sufficient statistics in R L and R R normalized by ∇ 2 ψ(β R ) + 2λΩ, further multiplied by the sample proportions in R L and R R . By selecting the candidate split that maximizes the quadratic term, we will end up with two subregions different in the sufficient statistics means and reasonably balanced in sample sizes.
To compute the quadratic approximation, we need subsample proportions n R s,L /n R , n R s,R /n R , average sufficient statisticsz R s,L ,z R s,R , and the inverse matrix of ∇ 2 ψ(β R )+2λΩ. For the candidate splits based on the same covariate, {n R s,L , n R s,R ,z R s,L ,z R s,R } can be computed efficiently for all split points by scanning through the samples in R once. For all candidate splits, this takes O(dn R k) operations in total. The matrix ∇ 2 ψ(β R ) is shared by all candidate splits and needs to be computed only once. The difficulty is that ∇ 2 ψ(β) is often unavailable in closed form. However, since ∇ 2 ψ(β R ) is the covariance matrix of the sufficient statistics z(y) if the responses y are generated from the model parameterized by β R , we apply Lindsey's method to estimate β R and compute the covariance matrix of the sufficient statistics based on the multinomial cell probabilities, which takes O(k 2 B). Appendix F shows the resulting covariance matrix approximates ∇ 2 ψ(β R ) with a fine discretization. The total time complexity of the above splitting procedure is summarized in the following Proposition 3.
Proposition 3 Assume that there are S candidate splits, k basis functions, d covariates, B discretization bins, and n R observations in the current region. Then the splitting step for LinCDE trees is of time complexity O(dn R k + k 2 B + k 3 + Sk 2 ).
According to Proposition 3, the computation time based on the quadratic approximation is significantly reduced compared to running Lindsey's methods in R s,L and R s,R for all candidate splits, which takesÕ(S(n R k + k 2 B + k 3 )).
Having found the best split s max , we partition R into two subregions R smax,L and R smax,R , and repeat the splitting procedure in the two subregions. Along the recursively partitioning, the response distribution's heterogeneity is reduced. The fitting and the splitting steps of LinCDE trees are summarized below, the complete algorithm is given in Algorithm 1, and implementation details are displayed in Section 10. Stopping criteria for LinCDE trees are discussed in Appendix D.
• Fitting (LinCDE tree). At a region R: 1. Count the number of observations {n R,b } in each bin.
2. Fit a Poisson regression model with the response variable {n R,b }, regressors z(y b ), the offset log(κ(y b )), and a weighted ridge penalty. 5 Denote the estimated coefficients byβ R .
• Splitting (LinCDE tree). At a region R: 1. Compute {n R s,L , n R s,R ,z R s,L ,z R s,R } for each candidate split, and approximate ∇ 2 ψ(β R ) by z(y)'s covariance matrix in R usingβ R .
2. For each candidate split s ∈ S, compute the quadratic approximation ∆ (R, s) by Proposition 2, and choose the split s max = arg max s∈S∆ (R, s).

Algorithm 1: LinCDE tree
Start at the full covariate space. 1. Apply Fitting (LinCDE tree) and obtain the natural parameter estimatorβ. 2. Apply Splitting (LinCDE tree) and obtain the optimal split s max . 3. Repeat steps 1 and 2 to the left and right children of s max until the stopping rule is satisfied, e.g., the maximal tree depth is reached. Output the natural parameter estimatorβ in each subregion. To conclude the section, we demonstrate the effectiveness of LinCDE trees in three toy examples. We generate 10 covariates randomly uniformly on [−1, 1]. The response follows with three different local densities f l (x), 1 ≤ l ≤ 3, varying in variance, number of modes, and skewness. The response distribution is determined by the first two covariates and independent of the rest. In Figure 2, we plot the average conditional density estimates in the three subregions. LinCDE trees are able to distinguish the densities differing in the above characteristic properties and produce good fits. We also compute the normalized importance score-the proportion of overall improvement in the split-criterion attributed to each splitting variable. In all settings, the first two covariates contribute over 99% importance. In other words, LinCDE trees focus on the first two influential covariates and avoid splitting at nuisance covariates.

LinCDE Boosting
Although LinCDE trees are useful as stand-alone tools, our ultimate goal is to use them as weak learners in a boosting paradigm. Standard tree boosting (Friedman, 2001) builds an additive model of shallow trees in a forward stagewise manner. Though a single shallow tree is high in bias, tree boosting manages to reduce the bias by successively making small modifications to the current estimate. 6 We proceed with the boosting idea and propose LinCDE boosting. Starting from a null estimate, we iteratively modify the current estimate by modifying the natural parameter functions via a LinCDE tree. In particular, at the t-th iteration, Section 5.1 gives details. We remark that the LinCDE tree modifier γ t (x) for boosting is an expanded version of that in Section 4: in previous LinCDE trees, all samples share the same carrying density κ(y), while in LinCDE trees for boosting, the carrying densities κ(y)e z(y) β t (x)−ψ(β t (x)) differ across units. We elaborate on LinCDE trees with heterogeneous carrying densities in Sections 5.2 and 5.3. Before discussing the details of LinCDE boosting, we compare LinCDE boosting and LinCDE trees on a toy example in Figure 3. We consider a locally Gaussian distribution with heterogeneous mean and variance 6. We remark that another successful ensemble method -random forests (Breiman, 2001)-are not appropriate for LinCDE trees. Random forests construct a large number of trees with low correlation and average the predictions. Deep trees are grown to ensure low-bias estimates, which is, however, unsatisfactory here because deep LinCDE trees will have leaves with too few observations for density estimation.  (14). The local densities are Gaussian with variances σ 2 ∈ {0.25, 1, 4} (first row), Gaussian mixtures with 1, 2, 3 components (second row), and Beta distributions with skewness κ ∈ {−0.8, 0, 0.8} (third row). In each trial, we sample 400 observations. We use 10 natural cubic splines , 5 degrees of freedom, and a maximal tree depth at 2. We repeat 100 times and plot the average fits against the true densities.
The covariate generating mechanism is the same as in Figure 2. We plot the average estimated conditional densities plus and minus one standard deviation. Though both LinCDE trees and LinCDE boosting produce good fits in all settings, the estimation bands of LinCDE  (16). We pick 4 landmarks corresponding to different conditional means and variances. From top to bottom, the conditional means decrease from 0.5 to −0.5. From left to right, the conditional variances increase from 0.6 to 1.6. In each trial, we sample 400 observations from the target distribution. We repeat each setting 100 times, and plot the average estimated conditional densities plus and minus one standard deviation against the true densities.
boosting are always narrower by around a half. The observation implies LinCDE boosting is more stable than LinCDE trees.

Additive Model in the Natural Parameter Scale
In LinCDE boosting, we build an additive model in the natural parameter scale of the density (11). We find a sequence of LinCDE tree-based learners with parameter functions {γ t (x)} 0≤t≤T −1 , and aggregate those "basis" functions to obtain the final estimate 7 In other words, at the t-th iteration, we tilt the current conditional density estimate 7. To stabilize the performance, we may shrink γ t (x) by some learning rate η ∈ (0, 1], and let β T (x) = is the updated normalizing function (depending on x). Appendix E shows that if the true conditional density is smooth, the approximation error of LinCDE boosting's function class (18) with splines z(y) will vanish as the number of splines k increases. We determine the LinCDE tree modifiers in (17) by log-likelihood maximization. We aim to find the modifier that produces the largest improvement in the objective Compared to the objective (12) of LinCDE trees, the only difference in (19) is the normalizing function φ β t (x) (γ t (x)). When β t (x) is a constant function, the normalizing functions of LinCDE trees and LinCDE boosting coincide. In Section 5.2 and 5.3, we demonstrate how LinCDE boosting's heterogeneous normalizing function complicates the fitting and splitting steps and propose corresponding solutions.

Fitting Step
For fitting, given a subregion R, the problem (19) can not be solved by Lindsey's method as in LinCDE trees, because β t (x) could be non-constant in the subregion R. Explicitly, instead of the single constraint in (5), we could have up to n R constraints As a result, the Lagrangian function (6) as well as subsequent discrete approximations for Lindsey's method are invalid. Fortunately, we can solve the fitting problem iteratively (Fitting (LinCDE boosting) below). Define bin probabilities We feed the marginal cell probabilitiesp b (R; β t (x)) to the fitting step as the baseline for modification. In Step 1, Lindsey's method produces a natural parameter modifier and a universal intercept for all samples in R. The intercept produced by Lindsey's method guarantees that the marginal cell probabilities to sum to unity, but not for every individual x i . In Step 2, we update the individual normalizing constants to ensure all constraints (20) are satisfied. In Proposition 4, we show that the fitting step of LinCDE boosting converges to the maximizer of the objective (19).
Proposition 4 Assume that λ = 0 and Y is supported on the midpoints {y b }, then the fitting step of LinCDE boosting converges, and the output γ t R satisfies We offer some intuition behind Proposition 4; i.e., why the fitting step of LinCDE boosting will stop at the likelihood maximizer. If β t (x) + γ is already optimal, then the average sufficient statistics z(y) under marginal probabilities (21) should match the observations, which yields the KKT condition of the Poisson regression in Lindsey's method. As a result, Lindsey's method will produce zero updates, and the algorithm converges.

Splitting Step
Reminiscent of the splitting step for LinCDE trees, we seek the split that produces the largest improvement in the objective (19). Proposition 2 is not valid due to the heterogeneity in β t (x), and we propose an expanded version.
Proposition 5 In a region R, let n R be the sample size and γ t R be the optimal update. Define the average sufficient statistics residuals as Given a candidate split s, define Then the improvement of the unpenalized conditional log-likelihood satisfies 2 ).
The average sufficient-statistic residualsr t R measures the deviation of the current estimator from the observations. By maximizing the quadratic approximation in Proposition 5, we find the candidate split s such thatr t R s,L andr t R s,R are far apart, and modify the current estimator differently in the left and right children determined by the selected split. The updated splitting procedure is summarized in Splitting (LinCDE boosting).
• Splitting (LinCDE boosting). In a region R: (22) by the average covariance matrix of z(y) in R. 2. For each candidate split s ∈ S, compute the quadratic approximation ∆ (R, s) by Proposition 5, and choose the split s max = arg max s∈S∆ (R, s).
The computation of the quadratic approximation is largely the same, except that the normalization matrix Ψ t s (γ t R ) varies across candidate splits and requires separate computation. To relieve the computational burden, we propose the following surrogate independent of candidate splits 8Ψ Proposition 6 gives the computational time complexity of the splitting procedure (Splitting (LinCDE boosting)). The computation time scales linearly with regard to the sample size multiplied by dimension and the number of candidate splits. The extra computation compared to LinCDE trees comes from residual calculations and individual normalizations.
Proposition 6 Assume that there are S candidate splits, then the splitting step for LinCDE boosting is of computational time complexityÕ(dn R kB + n R k 2 B + k 3 + Sk 2 ).

Pretreatment
In this section, we discuss two pretreatments: response transformation and centering. The pretreatments are helpful when the response is heavy-tailed and when the conditional distributions f y|x (y | x) vary wildly in location.

Response Transformation
Heavy-tailed response distributions are common in practice, such as income and waiting time. If the response is heavy-tailed, then in Lindsey's method, most bins will be approximately empty. As a result the model tends to be over-parameterized and the estimates tend to overfit.
8. In practice, we add a universal diagonal matrix to 1 to stabilize the matrix inversion. 9. The initialization β 0 (x) is usually a constant vector, e.g., zero vector, independent of the covariates.

Algorithm 2: LinCDE boosting
Initialize the natural parameter function β 0 (x). 9 for t = 1:T do 1. Apply Algorithm 1 with Fitting (LinCDE boosting), Splitting (LinCDE boosting), and obtain the optimal LinCDE tree modifierγ t−1 (x). In response to the heavy-tailed responses, we recommend transforming the response first to be marginally normally distributed. Often log and cube-root transformations are useful. In a more principled way, we can apply the Box-Cox transformation to the responses and choose the optimal power parameter to produce the best approximation of a Gaussian distribution curve. Once the model is fit to the transformed data, we map the estimated conditional densities of the transformed responses back to those of the original observations. See Figure 4 for an example.

Centering
For a distribution whose conditional components differ wildly in location, LinCDE needs a large number of sufficient statistics to capture local distributional characteristics. For instance, Figure 5 displays a conditional Gaussian mixture with location shift When we apply LinCDE boosting with k = 10 sufficient statistics, the estimates do not reproduce the bimodalities due to a lack of flexibility. We call this the "disjoint support" problem.
A straightforward solution to the disjoint support problem is to increase the number k so that the sufficient statistics z(y) are adequately expressive. As a consequence, the number of components in the parameter function β(x) goes up. This approach is prone to overfitting, especially when there are a small number (∼ 20) samples in a terminal node. In addition, this approach will significantly slow down the splitting procedure, which scales O(k 3 ) by Proposition 6.
Our solution is to center the response prior to fitting the LinCDE model. Since the difference in location causes the disjoint support problem, we suggest aligning the centers of the conditional densities in advance. Explicitly, we first estimate the locations via some conditional mean estimator and then subtract the estimates from the responses. The support of the residuals are less heterogeneous, and we apply LinCDE boosting to these residuals to capture additional distributional structures. Finally, we transform the resulting density estimates back to those of the responses. The procedure is summarized in Algorithm 3.
Algorithm 3: Centering 1. Estimate the conditional meanĥ(x) using the training data {(x i , y i )}. Compute the residuals r i = y i −ĥ(x i ).

Apply LinCDE boosting to
Centering splits the task of conditional distribution estimation into conditional mean estimation and distributional property estimation. For centering we have available a variety of popular conditional mean estimators, such as the standard random forest, boosting, and neural networks. Once the data are centered, LinCDE boosting has a more manageable task. Figure 5 shows that with centering, LinCDE boosting is able to reproduce the bimodal structure in the above example with the same set of sufficient statistics.

Simulation
In this section, we demonstrate the efficacy of LinCDE boosting on simulated examples. truth LinCDE boosting Figure 5: Conditional density estimation with and without centering. We consider the conditional density (23) and pick 3 landmarks corresponding to different locations. The first row plots LinCDE boosting's estimates without centering, and the second row plots the estimates augmented with true means. In each trial, we sample 1000 observations from the target distribution. We repeat each setting 100 times, and plot the average estimated conditional densities. In both settings, LinCDE boosting uses k = 10 sufficient statistics and 20 response bins.

Data and Methods
Consider d = 20 covariates randomly generated from uniform [−1, 1]. The responses given the covariates are sampled from the following distributions: • Locally Gaussian distribution (LGD): At a covariate configuration, the response is Gaussian with the mean determined by x (1) and x (2) , and the variance determined by x (2) . Covariates x (3) to x (20) are nuisances variables; • Locally Gaussian or Gaussian mixture distribution (LGGMD): where the means and variances are The mean is determined by x (1) . The modality depends on x (2) : in the subregion x (2) ≥ 0.2, the response follows a bimodal Gaussian mixture distribution, while in the complementary subregion, the response follows a unimodal Gaussian distribution. The skewness or symmetry is controlled by x (3) in the Gaussian mixture subregion: larger absolute values of x (3) imply higher asymmetry. Overall, the conditional distribution has location, shape, and symmetry dependent on the first three covariates. Covariates x (4) to x (20) are nuisance variables.
The training data set consists of 1000 i.i.d. samples. The performance is evaluated on an independent test data set of size 1000. We compare LinCDE boosting with quantile regression forest and distribution boosting. 10 There are a number of tuning parameters in LinCDE boosting. The primary parameter is the number of trees (iteration number). Secondary tuning parameters include the tree size, the learning rate, and the ridge penalty parameter. On a separate validation data set, we experimented with a grid of secondary parameters, each associated with a sequence of iteration numbers, and select the best-performing configuration. By default, we use k = 10 transformed natural cubic splines and a Gaussian carrying density We use a small learning rate η = 0.01 to avoid overfitting. We use 40 discretization bins for training, and 20 or 50 for testing. The simulation examples do not have heavy-tail or disjoint support issues, and thus no pretreatments are needed.

Results of Conditional Density Estimation
Let the oracle be provided with the true density, and the null method estimates a marginal Gaussian distribution. We consider the following metric where · denotes the test conditional log-likelihood of a specific method. The criterion is analogous to the goodness-of-fit measure R 2 of linear regression. It measures the performance of the method relative to the oracle; larger values indicate better fits, and the ideal value is one. Quantile regression forests and distribution boosting estimate conditional quantiles instead of densities. To convert the quantile estimates to density estimates, we define a grid of bins with endpoints y b,L and y b,R , and approximate the density in bin b bŷ whereq −1 (y) represents the inverse function of the quantile estimates. As the bin width shrinks,f b is less biased but of larger variance. In simulations, we display the results with 20 bins and 50 bins (Appendix G). We observe that LinCDE boosting is robust to the bin size, while distribution boosting and quantile regression forests prefer 20 bins due to the smaller variances. Figure 6 presents the goodness-of-fit measure (24) of the three methods under the LGD and LGGMD settings. In both settings, LinCDE boosting leads in performance, improving the null method by 60% to 80% of the oracle's improvements. changes, and in scale as x (2) changes; in the LGGMD setting, the estimated conditional densities vary in location as x (1) changes, in shape as x (2) changes, and in symmetry as x (3) changes. To further illustrate the ability of LinCDE boosting to detect influential covariates, we present the importance scores in Figure 9. In the LGD setting, LinCDE boosting puts around 87% of the importance on x (1) and x (2) , while quantile regression forest distributes more importance on the nuisances (x (1) and x (2) accounting for 40%). In the LGGMD setting, LinCDE boosting is able to detect all influential covariates x (1) , x (2) , x (3) , while the quantile regression forest only recognizes x (1) .

Results of Conditional CDF Estimation
Here we evaluate the conditional CDF estimates of the three methods, bringing the comparisons closer to the home court of distribution boosting. We consider the average absolute where {u j } is an evenly spaced grid on [0, 1], and q(u | x) denotes the u quantile at the covariate value x. To compute the CDF estimates, for distribution boosting and quantile regression forest, we directly invert the estimated quantiles to CDFs. For LinCDE boosting, we compute the multinomial cell probabilities with a fine grid (50 bins) and obtain the CDFs based on the cell probabilities. LGGMD (right panel) settings. We normalize the importance scores to sum to one. In the LGD setting, both methods detect all the influential covariates x (1) and x (2) . In the LGGMD setting, LinCDE boosting identifies all the influential covariates x (1) , x (2)) , and x (3) , while quantile regression forest only identifies x (1) . Figure 10 depicts the AAE metrics. In both settings, LinCDE boosting produces the smallest AAE. Notice that Though distribution boosting and quantile regression forest estimate the quantiles well, the CDF estimates can be harmed by the implicit density estimator multiplied. In Appendix G, we also compare the CDF estimates using Cramér-von Mises distance and observe consistent patterns to what we see here.

Results of Conditional Quantile Estimation
Here the comparisons are in the home court of quantile random forests. We evaluate the conditional quantile estimates of the three methods. We compute the quantile losses at {5%, 25%, 50%, 75%, 95%} levels (Table 1). For LinCDE boosting, we compute the multinomial cell probabilities (50 bins) and obtain the quantiles based on the cell probabilities. Despite the fact that quantile-based metrics should favor quantile-based methods, we observe that the performance of LinCDE boosting is similar.
We also construct 50% and 90% prediction intervals based on the quantiles. In Table 1, we display the coverages and widths of the 90% prediction intervals. All methods produce fairly good coverages. Results of 50% prediction intervals are consistent and can be found in Appendix G.

Computation Time
We compare the computation time of the three methods. 11 We normalize the training time by the number of trees. In Table 2

Real Data Analysis
In this section, we analyze real data sets with LinCDE boosting. The pipeline is as follows: first, we split the samples into training and test data sets; next, we perform 5-fold crossvalidation on the training data set to select the hyper-parameters; finally, we apply the estimators with the selected hyper-parameters and evaluate multiple criteria on the test data set. We repeat the procedure 20 times and average the results. As for the centering, we use random forests as the conditional mean learner.

Old Faithful Geyser Data
The Old Faithful Geyser data records the eruptions from the "Old Faithful" geyser in the Yellowstone National Park (Azzalini and Bowman, 1990) and represents continuous measurement from August 1 to August 15, 1985. The data consists of 299 observations and 2 variables: eruption time and waiting time for the eruption. We estimate the conditional distribution of the eruption time given the waiting time.
In Figure 11, we plot the eruption time versus the waiting time. There is a clear cutoff at 70min: for any waiting time over 70min, the distribution of eruption time is bimodal, while for any waiting time below 70min, the distribution is unimodal. In Figure 11, we display the estimated conditional densities of LinCDE boosting at waiting time 85min and 60min. LinCDE boosting is capable of detecting the difference in modality. In Table 3, we summarize the comparison between LinCDE boosting, distribution boosting, and quantile regression forest regarding the negative log-likelihoods and quantile losses. LinCDE boosting outperforms in log-likelihood, and is competitive in quantile losses. We remark the AAE and Cramér-von Mises distance can not be computed since we do not have the true conditional distributions.

Human Height Data
The human height data is taken from the NHANES data set: a series of health and nutrition surveys collected by the US National Center for Health Statistics (NCHS). We estimate the conditional distribution of the standing height. We consider two subsets: 542 samples in the age range 14 to 17, and 1956 samples in the age range 14 to 40. In the smaller subset, we only consider two covariates: age and poverty; in the larger subset, we consider 9 covariates, including age, poverty, race, gender, etc. In the smaller subset, we tune by cross-validation; in the larger subset, we split the data set for validation, training, and test (proportion 2 : 1 : 1), and tune on the hold-out validation data.
The distribution of heights combining male and female is used as a typical illustration of bimodality (Devore and Peck, 2005). However, Schilling et al. (2002) point out the separation between the heights of men and women is not large enough to produce the bimodality. In Figure 12, we demonstrate the histogram of heights of white teenagers in the age range 15-19. The distribution of the combined data is sightly bimodal. Boys' heights are larger and more concentrated. We also provide LinCDE boosting's conditional density estimates obtained from the larger data set. Overall, the estimates accord with the histograms. The estimates without the covariate gender is on the borderline of unimodal and bimodal. The estimates with the gender explains the formation of the quasi-bimodality.
The comparisons of log-likelihood and quantile losses are summarized in Table 3. In both the larger and the smaller data sets, LinCDE boosting performs the best concerning the log-likelihood. The advantage is more significant in the smaller data set. The reason is that in the larger data set, there are more covariates, and the conditional mean explains a larger proportion of the variation in response, while in the smaller data set, the conditional distribution after the centering contains more information, such as the bimodality, which can be learnt by LinCDE boosting.

Relative Spinal Bone Mineral Density Data
The relative spinal bone mineral density (spnbmd) data contains 485 observations on 261 North American adolescents. The response is the difference in spnbmd taken on two consecutive visits divided by the average. There are three covariates: sex, race, and age (the average age over the two visits). We estimate the conditional distributions of the spnbmd. The comparisons of log-likelihood and quantile losses are summarized in Table 3. LinCDE boosting performs the best concerning the log-likelihood. The scatterplot of spnbmd versus age demonstrates serious heteroscedasticity: the variances reach the climax at age 12 and decrease afterward. In Figure 13, we plot LinCDE boosting's estimates at age 12, 15 and 20 of white females. The spreads of the conditional densities decrease as the age grows. At age 12, the spnmd distribution is right-skewed, while those at age 15 and 20 are approximately symmetric.  Figure 12: LinCDE boosting applied to the height data. On the left, we plot the histogram of heights of white teenagers in the age range 15-19. On the right, we plot the average estimated density by LinCDE boosting from the larger data set at age 17, race white, and other covariates fixed at the corresponding medians. The estimated conditional density without the gender is on the borderline between unimodal and bimodal. Furthermore, we contrast the histograms of male and female heights against that of the combined data in the left plot, explaining the formation of the quasi-bimodality. We also depict the estimated densities of females, males (right panel), which accord with the histograms (left panel).

Discussion
In this paper, we propose LinCDE boosting for conditional density estimation. LinCDE boosting poses no specific parametric assumptions of the density family. The estimates reflect a variety of distributional characteristics. In the presence of unrelated nuisance covariates, LinCDE boosting is able to focus on the influential ones. So far, we have discussed only univariate responses. Multivariate responses emerge in multiple practical scenarios: locations on a 2D surface, joint distributions of health indices, to name a few. Lindsey's method and thus LinCDE boosting can be easily generalized to multivariate responses. Assume the responses are p-dimensional, multivariate LinCDE boosting considers the density (11) with sufficient statistics involving y (1) to y (p) . As an illustrative example, if we use sufficient statistics {y (i) } and {y (i) y (j) }, the resulting density will be a multivariate Gaussian (see Efron and Hastie, 2016, chap. 8.3 for the galaxy data example). The response discretization now divides the hyper-rectangle response support into equal-sized p-dimensional bins and the rest of LinCDE boosting procedures carry over. The cost of multivariate responses is the exponentially growing number of bins and sufficient statistics, which requires more samples as well as computational power. In contrast, conditional quantiles for multi-dimensional responses are relatively less straightforward but several promising proposals have been proposed (Barnett, 1976;Chaudhuri, 1996;Kong Figure 13: LinCDE boosting applied to the relative spinal bone mineral density data. On the left, we display the scatterplot of spnbmd versus age. There is serious heteroscedasticity. On the right, we plot LinCDE boosting's estimates at age 12, 15 and 20 of white females. The spreads of the conditional densities decrease as the age grows. At age 12, the spnbmd distribution is right-skewed, while those at age 15 and 20 are approximately symmetric. and Mizera, 2012;Carlier et al., 2016). We refer readers to Koenker (2017) and references therein.
There are several exciting applications of LinCDE boosting.
• Online learning. Online learning processes the data that become available in a sequential order, such as stock prices and online auctions. As opposed to batch learning techniques which generate the best predictor by learning on the entire training data set once, online learning updates the best predictor for future data at each step. Online updating of LinCDE boosting is simple: we input the previous conditional density estimates as offsets, and modify them to fit new data.
• Conditional density ratio estimation. A stream of work studies the density ratio model (DRM), particularly the semi-parametric DRM. The density ratio can be used for importance sampling, two-sample testing, outlier detection (see Sugiyama et al., 2012, for an extensive review). Boosting tilts the baseline estimate parametrically based on the smaller group.
One future research direction is adding adaptive ridge penalty in LinCDE boosting. Recall the fits in the LGGMD setting (Figure 8) where the conditional densities change from a relatively smooth Gaussian density to a curvy bimodal Gaussian mixture, the estimates get stuck in between: in the smooth Gaussian subregions, the estimates produce unnecessary curvatures; in the bumpy Gaussian mixture subregions, the estimates are not sufficiently wavy. The lack-of-fit can be attributed to the universal constraint on the degrees of freedom: the constraint may be too stringent in some subregions while lenient in others.

Software
Software for LinCDE is made available as an R package at https://github.com/ZijunGao/ LinCDE. The package can be installed from GitHub with 1 install . packages ( " devtools " ) 2 devtools :: install _ github ( " ZijunGao / LinCDE " , build _ vignettes = TRUE ) The package comes with a detailed vignette discussing hyper-parameter tuning and demonstrating a number of simulated and real data examples.

Appendix A. Regularization Parameter and Degrees of Freedom
We discuss how the hyper-parameter λ relates to the degrees of freedom. 12 According to Eq. (12.74) in Efron and Hastie (2016), the effective number of parameters in Poisson models takes the form where n b are responses andη b are estimates of the natural parameter. In the following Proposition 7, we obtain an approximation of Eq. (27) in the ridge Poisson regression via a quadratic expansion of the objective function.
Proposition 7 Assume the responses n b are generated independently from the Poisson model with conditional means µ * b = κ b e z b β * (z includes the intercept) and the corresponding natural parameter η * b = log(µ * b ). Let λ > 0 be the hyper-parameter of the ridge penalty. 13 • For arbitrary β, the second order Taylor expansion at β * of the Poisson log-likelihood with ridge penalty is for some constant C > 0 independent of β.
• Letβ be the minimizer of the quadratic approximation (28), then where H β * = Z W Z is the Hessian matrix of the negative Poisson log-likelihood evaluated at β * .
We prove Proposition 7 in Appendix F. As a corollary, if Ω is a scalar multiple of the identity matrix, Eq. (29) agrees with the degrees of freedom formula Eq. (7.34) in Hastie et al. (2009). In practice, β * is unknown and we plugβ in Eq. (29) to compute the number of effective parameters. We remark that Eq. (29) includes one degree of freedom for the intercept and in the LinCDE package we will use Eq. (29) minus one as the degrees of freedom.
12. The degrees of freedom are derived under the Poisson approximation (7) of the original likelihood. 13. λ = nλ for the λ in (8). We use 10 transformed natural cubic splines as basis functions and increase the degrees of freedom in Lindsey's method from 2 to 7 (left to right). Figure 14 plots the estimated log-densities of a Gaussian mixture distribution from Lindsey's method under different degrees of freedom. As the degrees of freedom (excluding the intercept) increase from 2 to 7, the log-density curves evolve from approximately quadratic to significantly bimodal. At 7 degrees of freedom, the estimated log-density is reasonably close to the underlying truth.

Appendix B. Linear Transformation for Basis Construction
We expand on the linear transformation used in the basis construction in Section 3. Let the eigen-decomposition of Ω be U DU , where U ∈ R k×k is orthonormal, and D ∈ R k×k is a diagonal matrix with non-negative diagonal values ordered decreasingly. Define the lineartransformed cubic spline basesz(y) = U z(y) and the corresponding coefficientsβ = U β. Thenz(y) β = z(y) β andΩ = U ΩU = U U DU U = D. Therefore, z(y) → U z(y) is the desired linear transformation.

Appendix C. Approximation Performance of Proposition 2
In Figures 15 and 16, we empirically demonstrate the efficacy of the quadratic approximation suggested by Proposition 2 . We let the conditional densities depend solely on x (1) and jump ( Figure 15) or vary continuously (Figure 16) in conditional variance, modality, or skewness. We observe that at candidate splits where the left child and the right child are similar, e.g., all the candidate splits based on x (2) , the quadratic form and the exact log-likelihood difference are almost the same. At candidate splits where the left child and the right child are different, e.g., x (1) = 0 in Figure 15, the exact log-likelihood difference is large, and the quadratic form is sufficiently close to the difference to determine the optimal split. We remark that to gain robustness, we set the quadratic approximation to zero if one of

log−likelihood difference quadratic approximation
Tree-structured conditional densities Figure 15: Comparison of the log-likelihood difference (13) and the quadratic approximation in Proposition 2. There are two subregions: x (1) ≤ 0 and x (1) > 0. In different subregions, the conditional densities are different in conditional variance (σ 2 ), modality, or skewness (each column corresponds to a type of difference). In the same subregion, the conditional density does not change. In each trial, we sample 100 observations and 5 covariates. We use 5 transformed natural cubic splines as basis functions, and 30 candidate splits for each covariate equally spread across their ranges. We plot the log-likelihood difference (red) and the quadratic approximation (black) for candidate splits of x (1) (influential) and x (2) (nuisance), respectively. The results are aggregated over 100 times.
the candidate split's children contain less than 10 samples, which leads to the imperfect approximation at the boundary candidate splits.

Appendix D. Stopping Criteria for LinCDE Trees
We discuss stopping criteria for LinCDE trees. There is no universally optimal choice. If we build a single tree learner, the preferred strategy, at least for regression and classification trees according to Breiman et al. (1984), is to grow a large tree, then prune the tree using Continuously varying conditional densities Figure 16: Comparison of the log-likelihood difference (13) and the quadratic approximation in Proposition 2. The conditional densities vary continuously with |x (1) | in variance (σ 2 ), modality, or skewness (each column corresponds to a type of difference). In each trial, we sample 100 observations and 5 covariates. We use 5 transformed natural cubic splines as basis functions, and 30 candidate splits for each covariate equally spread across their ranges. We plot the log-likelihood difference (red) and the quadratic approximation (black) for candidate splits of x (1) (influential) and x (2) (nuisance), respectively. The results are aggregated over 100 times.
the cost-complexity pruning. If we train a random forest learner, then Breiman (2001) recommends stopping the splitting process only when some minimum node size, default to be 5 in the package randomForest (Liaw and Wiener, 2002), is reached. As for tree boosting, Friedman (2001) uses trees with the same number of terminal nodes. The number of terminal nodes is treated as a hyper-parameter of the boosting algorithm and tuned to maximize the performance on the data set at hand. The aforementioned stopping criteria generalize to LinCDE trees straightforwardly. Two options are currently available in codes: (1) stop when the tree depth reaches some prefixed level; (2) stop when the decrease in the objective fails to surpass a certain number-a greedy top-down approach. We don't recommend the criterion of stopping until a terminal node is pure, because there will be insufficient samples at terminal nodes for density estimation.

Appendix E. Approximation Error for LinCDE Boosting
We characterize the expressiveness of LinCDE boosting's function class class (18) with splines z(y).

Proposition 8 Assume the log-conditional-density function
where C > 0 is a constant depending on B.
The proof of Proposition 8 can be found in Section F.

Appendix F. Proofs
In this section, we provide proofs of aforementioned propositions and claims.

F.1 Proof of Claim 1
Proof. For u ∈ R k , u lies in the null space of Ω z if and only if u Ωu = 0. By the definition of Ω, 0 = u Ωu = (z (y) u) 2 dy, which implies that (z(y) u) = z (y) u = 0 almost everywhere. On one hand, if z(y) u is linear or quadratic, then (z(y) u) is automatically zero everywhere. On the other hand, since z(y) are cubic spline bases, z(y) u is piece-wise cubic, and (z(y) u) = 0 implies that z(y) u is piece-wise linear or quadratic. Because z (y) and z (y) are both continuous, z(y) u is also second-order continuous, z(y) u must be linear or quadratic in y.

F.2 Proof of Proposition 7
Proof. For arbitrary µ b , η b such that η b = log(µ b ), the second order Taylor expansion at µ * b of the Poisson log-likelihood of the b-th sample (y b ; η b ) is By definition, for some constant C independent of η * b . Plug Eq. (32) into Eq. (31), where the remainder r(y (33), we sum over all samples and finish the proof of Eq. (28).
We different the quadratic approximation (28) with respect to β and obtain the score function We solve the score function to obtain Then plugη = K + Zβ into Eq. (27) and we get This gives Eq. (29). The Hessian of the negative Poisson log-likelihood at β * is

F.3 Proof of Proposition 2
Proof. We simplify the differences in the log-likelihood, where we use n R s,Lz R s,L + n R s,Rz R s,R = nz R . The score equation ofβ R implies and similarly forβ R s,L ,β R s,R . Plug Eq. (35) into Eq. (34), we obtain By the Taylor expansion of ψ(β) and ∇ψ(β) atβ R , Plug Eq. (37) into Eq. (36), we get Finally by Eq. (35) and the assumption that ∇ 2 ψ(β R ) + 2λΩ is invertible, Thenβ R s,L −β R z R s,L −z R and similarly for the right child. Plug Eq. (39) into Eq. (38), and we finish the proof.

F.4 Proof of Claim 9
Claim 9 Assume that in R, y is supported on the midpoints {y b }, then the covariance matrix approximation by Lindsey's method equals ∇ 2 ψ(β R ).
Proof. If y is supported on the midpoints {y b }, then the discretization in Lindsey's method is accurate. As a result, the estimator of Lindsey's method is the exact loglikelihood maximizerβ R . Furthermore, the multinomial cell probabilities based on the Linsey's method's estimator is indeed the response distribution indexed byβ R . Next, by Lehmann and Romano (2006), ∇ 2 ψ(β R ) equals the population covariance matrix of the sufficient statistics generated from the distribution indexed byβ R . Therefore, the covariance approximation by Lindsey's method equals ∇ 2 ψ(β R ).

F.5 Proof of Proposition 3
Proof. We order the observations and count responses in each bin (Õ(n R )). Next, we run Newton-Rapson algorithm. In each iteration, the computation of the gradient vector To compute the quadratic approximation, {n R s,L , n R s,R ,z R s,L ,z R s,R } can be computed in O(dn R k) by scanning through the samples once per coordinate. (If observations are not ordered beforehand, we will haveÕ(dn R k), whereÕ denotes the order up to log terms.) Adding the diagonal matrix Ω is cheap, and the matrix inversion takes O(k 3 ) operations. For a candidate split, the quadratic term in Proposition 2 further takes O(k 2 ) time. Choosing the optimal split takes O(S). In summary, the complexity is O(dn R k + k 2 B + k 3 + Sk 2 ).

F.6 Proof of Proposition 8
The results are built on the expressiveness of logspline density models and tree boosting models. We first develop an intermediate property for unconditional densities, and then combine it with the representability of tree boosting models with arbitrary depth and number of trees. We will use C to denote constants depending on B and the concrete values may vary from line to line.
For unconditional densities, let the true log-density be h(y) ∈ S s 2 and h (s) (·) 2 ≤ B. By Barron and Sheu (1991), there exists a function h k ∈ Ω k,ζ such that where we use δ = 1/(k − ζ + 2) ≤ 2/k. Next, we turn to conditional densities. For |h(y | x)| ∈ S s 2 and h(· | x) ∞ ≤ B, by Barron and Sheu (1991, Lemma 3), there exists a unique vector solution β(x) to the minimization problem We next show β(x) is measurable on [0, 1] d . In fact, define the mapping

F.7 Proof of Proposition 4
Proof. Without loss of generality, assume R is the covariate space. For simplicity, we omit the superscript of β t (x) and denote the natural parameter functions by β(x).
We first show γ * = arg max γ (R; β(x) + γ). For λ = 0, if ∆γ = 0, the KKT condition of the last Poisson regression is If in R, Y is supported on the midpoints {y b }, Therefore, by Eq. (42) and Eq. (43), which implies the KKT condition of maximizing (R; β(X) + γ) is satisfied at γ * . Since the log-likelihood (R; β(X) + γ) is strictly concave, then γ * is the unique maximizer. We next prove the algorithm converges. Plug in the MLE estimate of the intercept of the Poisson regression into the log-likelihood and we have By Jensen's inequality, Notice that Therefore, by Eq. (45), Eq. (46) and Eq. (47), By the strong concavity of poisson , there exists a universal δ > 0 such that Thus, Since the conditional log-likelihood (R; β(x) + γ) is bounded, thus can not increase linearly forever. Therefore, ∆γ 2 → 0, i.e. the algorithm converges.

F.8 Proof of Proposition 5
Proof. The proof follows that of Proposition 2. We simplify the difference of log-likelihoods, The score equation of γ t R implies and similarly for γ t R s,L , γ t R s,R . Plug Eq. (49) into Eq. (48), we obtain By the Taylor expansion of ψ(β) and ∇ψ(β), Plug Eq. (51) into Eq. (50), we get By Eq. (49), wherez t R s,L := x i ∈R s,L ∇ψ(β t (X i ) + λ t R )/n R s,L , and similarly forz t R s,R . Finally, by the assumption of invertibility, Plug Eq. (53) into Eq. (52), and we finish the proof.
F.9 Proof of Proposition 6 Proof. We first evaluate the complexity of the fitting step of LinCDE boosting. The computation of offset is O(n R kB), and we store the cell probabilities p b (β t (x i )). The first step of the fitting runs a penalized Lindsey's method and takes O(n R + k 2 B + k 3 ) as shown in the proof of 3. The second step in the fitting takes O(n R B + kB), and we update the cell probabilities to p b (β t (x i ) + γ t R ). It takes O(n R k 2 B) to compute the surrogate normalization matrixΨ t (γ t R ) and an extra O(k 3 ) for matrix inversion. It takesÕ(dn R kB) to compute all average residualsr t R . Finally, the quadratic approximation for all candidate splits and choosing the optimal one takes O(Sk 2 ). In summary, the time complexity isÕ(dn R kB + n R k 2 B + k 3 + Sk 2 ).

Appendix G. Additional Figures
In this section, we provide additional numerical results. In the left panel, the true density is a Gaussian mixture, and in the right panel, the true density is a beta distribution. In each trial, we sample 400 observations. We generate 10 natural cubic splines with knots equally spread across the range of observations, and tune the penalty parameter λ to achieve an effective 5 degrees of freedom. We repeat each setting 100 times and plot the average fits against the true densities. Figure 17 displays the performance of Lindsey's method in two toy examples. One target density is bimodal, and the other is skewed. In Lindsey's method, we use natural cubic splines, which are arguably the most commonly used splines because they provide good and seamless fits and are easy to implement, and transform them as discussed.

G.1 Additional Simulations
In both examples, the estimated densities of Lindsey's method match the true densities quite closely, except for tiny gaps at boundaries and peaks. Figure 18 presents the goodness-of-fit measure (24) of the three methods under the LGD and LGGMD settings with 50 bins. LinCDE boosting benefits from finer grids. In contrast, distribution boosting and quantile regression forest are hurt by larger numbers of bins due to higher variances, especially in the LGGMD setting where the densities themselves are bumpy.
Similar to AAE, another commonly used CDF-based metric is the cram/'er-von Mises distance (Friedman, 2019): where {u j } is an evenly spaced grid on [0, 1], and q(u | x) denotes the u quantile at the covariate value x. Figure 19 depicts the cram/'er-von Mises distance. In both settings, LGGMD (right panel) settings.
In Figure 20, we plot the coverages of the 50% prediction intervals. The results are consistent with those of 90% prediction intervals.

G.2 Air Pollution Data
The air pollution data (Wu and Dominici, 2020) focuses on PM2.5 exposures in the United States. 14 The responses are 3108 county-level PM2.5 exposures averaged from 2000 to 2018. We incorporate 16 weather, socio-economic, and demographic covariates, such as winter relative humidity, house value, and proportion of white people. The target is estimating the conditional density of the average PM2.5 exposure. We split the data into training, validation, and test (proportion 2:1:1), and tune on the hold-out validation data. We also identify influential covariates which may help find the culprits of air pollution and manage regional air quality.
The PM2.5 exposure varies vastly across counties. For example, the average PM2.5 of west coast counties may soar up to 12µg/m 3 due to frequent wildfires, while those of rural areas in Central America are typically below 8µg/m 3 . The difference in PM2.5 levels induces the disjoint support issue for LinCDE in Section 6, and thus we employ the centering. The comparisons of log-likelihood and quantile losses are summarized in Table 4. Centered LinCDE performs the best in log-likelihood and is comparable in quantile losses.
In Figure 21, we display how the estimated conditional densities change with respect to winter relative humidity and house value-top influential covariates identified by LinCDE.
• Winter relative humidity affects the locations of the conditional densities: as the humidity goes up, the PM2.5 concentration first increases, then decreases. One hy-pothesis of the inverted U-shaped relationship is that in dry counties, wind speed that is inversely correlated with the humidity is a powerful factor for PM2.5; in wet counties, moisture particles that accelerate the deposition process of PM2.5 are the driving force.
• House value is impactful to the scales of the conditional densities: more expensive houses are associated with more variable PM2.5 exposures. We conjecture that rural areas are generally low in PM2.5 while urban areas vary. Higher house values suggest a higher proportion of urban areas and thus less homogeneous air quality.    Figure 21: LinCDE boosting applied to the air pollution data. We plot the estimated densities of LinCDE boosting across different levels of top influential covariates: winter relative humidity (left panel) and house value (right panel). Other covariates are fixed at the corresponding medians. Winter relative humidity is influential to the locations of the conditional densities: as the humidity goes up, the PM2.5 concentration first increases, then decreases. House value is impactful to the scales of the conditional densities: more expensive houses are associated with more variable PM2.5 exposures.