Fit a function f on dataset X such that f(X) fits a histogram

I have dataset $X=\{\boldsymbol{x_1},\boldsymbol{x_2},\dots,\boldsymbol{x_n}\}$ and $Y=\{y_1,y_2,\dots,y_n\}$ and want to learn a function $f$ such that $y = f(\boldsymbol{x})$ can be approximated as much as possible (by whatever cost measure). The hard part is that due to the limitation of data source, I have $y_i$ only accessible as an interval $(y_i^{(1)}, y_i^{(2)})$ which covers the possible range of the real $y_n$, rather than an exact point. What's more, those intervals don't overlap and not necessarily have equal size. To be more specific, they looks like this:

(1001, 2000), (2001, 4000), (1001, 2000), (4001, 6000), ...

Now I wonder if there is any technique suitable for such situation. Can I learn such a function $f$ that the $f(X)$ coincides $Y$ as much as possible. For example, $P[f(\boldsymbol{x}\in (1001, 2000)]$ should resembles $P[y_n\in (1001, 2000)]$.

Hope I make myself clear. If not, please point out the obscureness.

• @user777 It seems KDE can only smooth the histogram. What I want is to is to learn a function from X to Y. It's not the same target. Feb 22 '16 at 16:28
• @user777 Smoothing the histogram is acceptable if it can help learn the mapping. But merely smoothing can't produce a mapping, I suppose. Feb 22 '16 at 16:35
• @user777 $x$ here is a vector rather than a scalar. Can KDE be applied to this case? I'm not familiar with KDE. Feb 22 '16 at 16:39

It seems that you want to do is to train a Gaussian process (GP) for regression with observation-dependent noise so as to approximate $f(\textbf{x})$. The idea is that your $y_i$ are noisy (uncertain) observations of the underlying function $f(\textbf{x})$, which has an unknown non-parametric form.
The standard observation noise model for GP regression is Gaussian. You did not specify how the intervals are obtained, so it's unclear if they are to be considered as confidence intervals of some kind or as hard bounds. Either way, as a fair starting point I would use your intervals to set the observation as $y_i = \frac{1}{2} \left(y_i^{(1)} + y_i^{(2)}\right)$ with associated observation noise as $\sigma_i = \alpha \left(y_i^{(2)} - y_i^{(1)}\right)$, where you set $\alpha$ depending on the interpretation of your intervals (e.g., I would set $\alpha$ between $\frac{1}{4}$ and $\frac{1}{2}$).