Non parametric estimation/regression for conditional distribution Context : one continuous variable $Y$ dependent on $X$ ($X$ can be anything)
Linear regression, generalized linear model... focus on estimating the conditional expectation $E(Y\mid X)$. I want to focus on estimating the conditional distribution.
Linear regression can be seen as a parametric model for the conditional distribution: $Y\mid X\hookrightarrow \mathcal{N}(X\beta,\sigma^2)$ with parameter $\beta$. It assumes that given $X$, $Y$ has a normal distribution.
I'm looking for methods where the distribution of $Y$ given $X$ is not assumed to belong to a predefined family of distributions : just having a certain regularity and maybe vaguely a certain shape. How this distribution is assumed to change when $X$ changes is something that may be predefined or not. Of course this may require quite a lot of training data. 
For non conditional distributions (just the distribution of $Y$, there is no $X$), I know a simple method : kernel density estimator.
What methods exist for the conditional case ?
 A: In a time series context you can take a look at non-parametric predictive distribution estimation for $p(X_{t+1} \mid X_t, X_{t-1}, \ldots)$. For this I suggest to look at hard LICORS, mixed LICORS, LICORS cabinet, and Hilbert space embeddings of predictive state representations -- to name a couple (disclaimer: I have been working on the first 2.)
As a pure regression framework, I just finished work on "Predictive State Smoothing (PRESS): Scalable non-parametric regression for high-dimensional data with variable selection", which makes almost no assumptions on the distributions for $y \mid \mathbf{X}$ and estimates an optimal Kernel smoother.  As others have already pointed to npreg or other kernel smoother implementations, I want to highlight that PRESS does not suffer from curse of dimensionality.  In the paper I show examples of a kernel smoother for an $\mathbf{X} \in R^{60,000 \times 784}$, i.e., $784$ covariates with $N = 60,000$ observations.  PRESS can estimate the kernel smoothing matrix for that within a couple of seconds. Traditional kernel smoothers break down (or just won't work) for more than a couple (5-ish) dimensions.
