Conditional Distribution Learning · Tutorial 623 of 989

Series Estimation of Conditional Densities

Build conditional density estimates from basis expansions, understand their approximation and estimation error, and check when the fitted function is a valid density.

Advanced 12 min read

What You'll Learn

Represent a conditional density using basis functions in the covariate and response.
Fit series coefficients from conditional-moment observations.
Separate truncation error from coefficient estimation error.
Verify normalization and identify why a series estimate may be negative.
Apply a simple positivity-preserving shrinkage to a fitted expansion.

Why Use a Series for a Conditional Density?

A conditional density describes how the distribution of a response changes with a covariate. Kernel methods estimate it locally; a series method instead builds it from a collection of basis functions. The basis gives a finite-dimensional representation, so estimation becomes a problem of fitting coefficients. The central design challenge is to make the representation flexible enough to capture changes in both the response and the covariate, while respecting the requirements of a density.

For a concrete setup, take $X$ in a covariate space with distribution $P_X$, and $Y$ in $[0,1]$. Suppose the conditional law of $Y$ given $X=x$ has density $f(y\mid x)$ with respect to Lebesgue measure. Assume this density is square-integrable under the product measure $P_X(dx)\,dy$. We use orthonormal basis functions $\{\phi_j\}_{j\geq 0}$ for $L^2(P_X)$ and $\{\psi_k\}_{k\geq 0}$ for $L^2([0,1])$, with $\phi_0(x)=1$ and $\psi_0(y)=1$. In particular, $\int_0^1\psi_k(y)\,dy=0$ for $k\geq1$.

Definition: A tensor-product series approximation to a conditional density is a finite sum of products of covariate and response basis functions. With the constant response term fixed, its normalized form is $$ f_{J,K}(y\mid x)=1+\sum_{j=0}^{J}\sum_{k=1}^{K}a_{jk}\phi_j(x)\psi_k(y). $$ The coefficients $a_{jk}$ describe how the conditional response-basis moments vary across the covariate. The indices $J$ and $K$ control the number of covariate and response basis functions retained.

The constant term is fixed at one because a conditional density must integrate to one for every $x$. If the response support is a different interval, the constant basis function must be scaled to integrate to one; the same principle applies. More generally, one can start with any response basis and enforce the normalization constraint on its coefficients. The form above builds that constraint in from the outset.

Fitting the Coefficients

For $k\geq1$, define the conditional response-basis moment $m_k(x)=E[\psi_k(Y)\mid X=x]$. Orthonormality gives the expansion $f(y\mid x)=1+\sum_{k\geq1}m_k(x)\psi_k(y)$ in squared-integrable terms. The coefficient of $\phi_j$ in $m_k$ is therefore $a_{jk}=E[\phi_j(X)\psi_k(Y)]$. Given independent observations $(X_i,Y_i)$, a direct coefficient estimate is the sample average $\hat a_{jk}=n^{-1}\sum_{i=1}^n\phi_j(X_i)\psi_k(Y_i)$.

This moment estimate is unbiased when the observations follow the stated joint distribution. It is especially simple when the basis functions are orthonormal under the true covariate distribution. In practice that distribution is usually unknown; one can instead fit regressions of $\psi_k(Y_i)$ on the covariate basis, using the sample Gram matrix. That version accounts for non-orthogonality in the observed design. Either way, the resulting conditional density estimate has the form $$ \hat f_{J,K}(y\mid x)=1+\sum_{j=0}^{J}\sum_{k=1}^{K}\hat a_{jk}\phi_j(x)\psi_k(y). $$

Worked Example: Estimating a Constant Conditional Moment

Let $\psi_1(y)=\sqrt{3}(2y-1)$, which has integral zero and squared integral one on $[0,1]$, and use only $\phi_0(x)=1$. For four observed responses $0.25, 0.75, 0.50, 0.25$, calculate the estimated coefficient.

The basis values are $-\sqrt{3}/2,\ \sqrt{3}/2,\ 0,\ -\sqrt{3}/2$. Thus $$ \hat a_{01}=\frac{1}{4}\left(-\frac{\sqrt{3}}{2}+\frac{\sqrt{3}}{2}+0-\frac{\sqrt{3}}{2}\right) =-\frac{\sqrt{3}}{8}\approx -0.2165. $$

The fitted function is $1-(\sqrt{3}/8)\psi_1(y)$, independent of $x$. Since $\psi_1$ ranges from $-\sqrt{3}$ to $\sqrt{3}$, its minimum fitted value is $1-3/8=0.625$. It is nonnegative and integrates to one, so this particular fitted function is a valid density for every $x$.

Normalization, Positivity, and Approximation

The fixed constant term ensures normalization, but positivity is a separate constraint. For arbitrary fitted coefficients, the sum of basis functions can dip below zero at some response values. Thus a normalized series is not necessarily a conditional density. Checking this distinction is essential before interpreting the fitted function as a probability model.

Theorem: For any finite real coefficients $a_{jk}$, the series $f_{J,K}(y\mid x)$ integrates to one over $y\in[0,1]$ for every $x$. It is a conditional density only if it is also nonnegative for almost every $y$, for the covariate values at which it is used.

Proof. Integrate the finite sum term by term. For each $k\geq1$, $\int_0^1\psi_k(y)\,dy=0$, so every nonconstant response term vanishes: $$ \int_0^1 f_{J,K}(y\mid x)\,dy =1+\sum_{j=0}^{J}\sum_{k=1}^{K}a_{jk}\phi_j(x)\int_0^1\psi_k(y)\,dy =1. $$ The sum is finite, so no limit or interchange issue arises. Nonnegativity is not implied by this calculation and must be imposed or checked separately. $\square$

Worked Example: A Normalized Series That Is Not a Density

Take $\hat f(y\mid x)=1+0.7\psi_1(y)$, with the same $\psi_1(y)=\sqrt{3}(2y-1)$. Its integral is one because $\psi_1$ integrates to zero. But at $y=0$, $\psi_1(0)=-\sqrt{3}$, so $$ \hat f(0\mid x)=1-0.7\sqrt{3}\approx -0.2124. $$

The fitted function is therefore not a density despite being normalized. One simple repair is to shrink the nonconstant part. If $g(y\mid x)$ has zero integral in $y$, then $1+\lambda g(y\mid x)$ remains normalized. When $M=\sup_{x,y}|g(y\mid x)|$ is finite and positive, taking $\lambda=1/(2M)$ ensures $1+\lambda g\geq1/2$. Here $M=0.7\sqrt{3}$, so $\lambda=1/(1.4\sqrt{3})\approx0.4124$; the shrunk coefficient is about $0.2887$, and the minimum becomes $0.5$. This guarantees validity but changes the fitted coefficients and is not generally the best approximation.

A principled alternative is constrained coefficient fitting: minimize a chosen data-fit criterion subject to $1+\sum_{jk}a_{jk}\phi_j(x)\psi_k(y)\geq0$ throughout the target domain, while retaining the zero-integral response basis. These are linear inequalities in the coefficients at each fixed $(x,y)$. Enforcing them everywhere can be difficult; checking only a finite grid may miss negative values between grid points. Other constructions use transformations or positive basis functions, but they change the fitting problem.

Theorem: Suppose the bases are orthonormal under $P_X$ and Lebesgue measure on $[0,1]$, respectively. Suppose $|\phi_j(x)|\leq A$ and $|\psi_k(y)|\leq B$ for all retained indices and almost every input. Estimate each coefficient by the sample average above, and let $f_{J,K}$ be the orthogonal projection of $f$ onto the retained tensor-product terms together with the constant function. Then $$ E\left[\|\hat f_{J,K}-f\|_{L^2(P_X\times dy)}^2\right] =\|f_{J,K}-f\|_{L^2(P_X\times dy)}^2+ \sum_{j=0}^{J}\sum_{k=1}^{K}\operatorname{Var}(\hat a_{jk}) \leq \|f_{J,K}-f\|_{L^2(P_X\times dy)}^2+\frac{(J+1)K A^2B^2}{n}. $$

Proof. Orthonormality makes the retained products $\phi_j\psi_k$ mutually orthonormal, and they are orthogonal to the constant function for $k\geq1$. The difference $f_{J,K}-f$ is orthogonal to every retained product by the projection property. The estimation error $\hat f_{J,K}-f_{J,K}$ lies in their span. Consequently, for each sample, $$ \|\hat f_{J,K}-f\|^2=\|\hat f_{J,K}-f_{J,K}\|^2+\|f_{J,K}-f\|^2. $$ The expected squared estimation error is the sum of the coefficient variances, since the squared norm in an orthonormal basis is the sum of squared coefficients. For a fixed pair $(j,k)$, the estimated coefficient averages $Z_i=\phi_j(X_i)\psi_k(Y_i)$, so independence gives $$ \operatorname{Var}(\hat a_{jk})=\frac{\operatorname{Var}(Z_1)}{n} \leq\frac{E[Z_1^2]}{n}\leq\frac{A^2B^2}{n}. $$ There are $(J+1)K$ retained coefficients. Summing this bound proves the inequality. The first term is truncation or approximation error; the second is coefficient estimation error. $\square$

Worked Example: A Covariate Term Omitted by Truncation

Let $X$ be uniform on $[0,1]$, set $\phi_1(x)=\sqrt{12}(x-\tfrac12)$, and use $\psi_1(y)=\sqrt{3}(2y-1)$. Consider $$ f(y\mid x)=1+0.3\phi_1(x)\psi_1(y). $$

Both basis functions have squared integral one and mean zero on their respective intervals. Their absolute maxima are $\sqrt{3}$, so $|0.3\phi_1(x)\psi_1(y)|\leq0.3(3)=0.9$. Hence $f\geq0.1$, and the zero-integral response term ensures it integrates to one. If we retain $K=1$ but set $J=0$, the only possible coefficient is $a_{01}=E[\phi_0(X)\psi_1(Y)]$. Under the displayed model this is $0.3E[\phi_1(X)]=0$. The resulting projection is the uniform density $1$, which misses the covariate-dependent term.

The integrated squared approximation error is $$ \int_0^1\int_0^1\left(0.3\phi_1(x)\psi_1(y)\right)^2\,dy\,dx =0.09\left(\int_0^1\phi_1(x)^2dx\right)\left(\int_0^1\psi_1(y)^2dy\right) =0.09. $$ Retaining $J=1,K=1$ represents the density exactly, so this truncation error becomes zero. This example isolates approximation error: increasing the basis dimensions improves representational capacity, while estimating the added coefficients introduces additional sampling variability.

Choosing and Interpreting a Series Fit

The error decomposition clarifies the role of the truncation dimensions. Small $J$ or $K$ can omit meaningful variation, while large values require estimating more coefficients. The theorem gives a transparent finite-sample bound under bounded bases and known orthonormality; it does not by itself prescribe optimal dimensions or guarantee a density. Those decisions depend on smoothness, sample size, the covariate design, and the loss used for evaluation.

In applications, inspect the fitted function for negative values over the region where predictions will be used, and verify its integral numerically as a diagnostic. Numerical checks do not prove nonnegativity everywhere. If positivity is essential, use an explicitly constrained or positivity-preserving procedure, then assess how that constraint changes the fit. Also distinguish integrated approximation quality from pointwise behavior: a small integrated squared error need not prevent a localized negative dip.

Common Pitfall

A frequent error is to treat “integrates to one” as synonymous with “is a density.” The zero-integral response basis guarantees normalization, but coefficients can still combine to produce negative values. The correction is to verify both conditions: unit integral in the response for each relevant covariate, and nonnegativity almost everywhere. Another error is to interpret increasing the basis size as an unqualified improvement. It reduces truncation error only when the larger span adds useful directions; it also adds coefficients to estimate and can increase variability.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

Why does a series using only response basis functions with zero integral preserve conditional normalization?
What additional condition must a normalized series satisfy to represent a conditional density?
Under the theorem’s assumptions, how does the coefficient-variance bound change if the number of retained response terms doubles while all else stays fixed?
In the truncation example, why is the $J=0$ projection uniform even though the true conditional density depends on $x$?
Why can checking nonnegativity on a finite grid fail to establish nonnegativity over the full covariate-response domain?

Back to Tutorials