# Regression of family of marginal density functions

From a large sample of triples $(X, Y, U)$ I need to estimate a function $(x, y, u) \mapsto f(x, y, u)$, such that for each fixed $x, y$, the function $u \mapsto f(x, y, u)$ is a density function; namely, the marginal density function of $U$, given $X = x$ and $Y=y$.

Is there a standard way of doing this? My initial thought was to do something along the lines of

1) For a collection of points in the $(x, y)$-grid, estimate (e.g. by some standard kernel density technique) $u \mapsto \widehat{f}_{x, y}(u)$ based on data points $(\tilde{x}, \tilde{y}, \tilde{u})$ where $\tilde{x}, \tilde{y}$ resides in some suitable neighbourhood about $(x, y)$ and $\tilde{u}$ varies freely.

2) Bilinearly interpolate between these points to obtain $f$.

Is this a reasonable approach?

It's also implemented e.g. in npcdens in the R np package or KDEMultivariateConditional in the Python statsmodels package.