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.
 A: 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}$). 
As a practical aside, note that many GP packages do not support observation-dependent noise (instead, they assume that the observation noise is the same for all observations). You might have to modify the code of the GP inference in order to do that, or you might want to code your own GP from scratch (it's not that difficult and it's very instructive).
As a second step, if you find the assumption of Gaussian observation too strong (either based on your fits, or because of prior knowledge of how the intervals are measured in your data), you might assume instead a uniform (boxcar) likelihood, or some other more appropriate likelihood. In this case inference in the Gaussian process is not exact any more, and you'll have to use approximations such as MCMC, EP (expectation-propagation), etc.
If you are not familiar with Gaussian Processes, the book by Rasmussen and Williams (2006) is the bible.

See also this answer for an alternative solution to a similar problem, in which however the lower/upper bounds are not jointly specified for all observations (some observations only have lower or upper bounds).
