Curve smoothing in the presence of non-gaussian uncertainty What options are available for smoothing 2-dimensional real data for which the the ordinate points are real intervals of the form
$(x_j , [y_{j0} , y_{j1}])$
In my case, the data is vague because of intrinsic measurement precision, and I want the knowledge of the quantity to be characterized by a uniform distribution over the interval $[y_{j0} , y_{j1}]$ (rather than Gaussian).
One option is to use a smoothing spline where the ordinate data is taken to be the means of the intervals and the points are weighted inverse-proportionally to the the size of the interval (in the Gaussian case, I guess this would be standard deviation instead). Unfortunately, this encourages the curve to pass through the middle of the intervals despite the fact that I have no reason to favour this.
What other options are available?
 A: A possible formulation is through gaussian processes. The unknown
smooth function $\eta(x)$ is seen as a stochastic process $Y(x)$ with
a given distribution that can be seen as a functional prior for
$\eta(x)$.  Then the estimation $\widehat{\eta}(x)$ can be the
posterior mean, i.e. the expectation of $Y(x)$ conditional on the set
of inequalities $y_{i0} \le Y(x_i)$ and $Y(x_i) \le y_{i1}$. However,
it is not evident to find a software computing this efficiently.
For a more practical solution,  the mid-point $y_i:=[y_{i0}+y_{i1}]/2$ can be taken as a 
response at $x=x_i$ as you did. Then most least-squares based smoothing 
methods can be adapted
by solving a Quadratic Programming (QP) problem instead of a 
least squares problem: since QP routines
are widely available, a program  will not be difficult to write. 
For instance, a constrained smoothing spline
can be found by QP. Let $\mathbf{y}$ be the vector of the $n$ mid-points $y_i$,
 and $\boldsymbol{\eta}$ be the vector of the unknown $\eta_i:=\eta(x_i)$.
In the usual spline smoothing, the estimate $\widehat{\boldsymbol{\eta}}$ is found by 
the minimization
$$
 \min_{\boldsymbol{\eta}} \: p \,\|\mathbf{y}-\boldsymbol{\eta}\|^2 
  + (1-p) \,\boldsymbol{\eta}^{\mathrm{T}} \mathbf{M} \boldsymbol{\eta}
$$
where $\mathbf{M}$ is a $n\times n$ matrix with rank $n-2$ depending on the design
points $x_i$, and $0 <p < 1$
is a smoothing parameter.
For the constrained spline, we add the two sets of $n$ constraints:
 $\boldsymbol{\eta} \ge \mathbf{y}_0$ and $\boldsymbol{\eta} \le \mathbf{y}_{1}$.
The unknown vector $\boldsymbol{\gamma}$ 
of $n-2$ "coefficients" i.e. of second order 
derivatives at interior nodes, which is needed e.g. to interpolate 
is related to $\boldsymbol{\eta}$ through $\boldsymbol{\eta} = \mathbf{K}\boldsymbol{\gamma}$ where the matrix $\mathbf{K}$
with dimension $n \times (n-2)$ is found (as well as $\mathbf{M}$) in the literature on smoothing splines. 
At least in a first approach, the parameter $p$ can be guessed.
When $p \approx 0$, the smoothed curve will be allowed to depart from the 
mid-points and  to be close to either end-points to reach a greater 
level of smoothing. However $p>0$ is needed to have a positive definite matrix
in the QP.  The quadprog R package can be used for, say, $n \le 100$,
and larger problems can probably be decomposed in blocks.
A: Kernal smoothing is possible but the probably special kernels should be used to limit the neighborhood of smoothing since you know that the noise term has a limited range.  
A: 
I could also have another source of noise which means that I could have the true curve appearing outside the interval.

In this case I suggest modeling the error by a generalized normal distribution, whose pdf is given by
$f_X(x; \alpha, \beta, \mu) \equiv \dfrac{\beta}{2\alpha\Gamma(1/\beta)} \; e^{-(|x-\mu|/\alpha)^\beta}$ 
where $\Gamma$ denotes the gamma function, and looks like

As you can see, for large values of the shape parameter $\beta$, the density is flat around the mean (zero in this example), approximating a uniform distribution. (Indeed, it converges pointwise to a uniform distribution inside the $(\mu-\alpha,\mu+\alpha)$ as $\beta \to \infty$.) Moreover, it is nonzero outside $(\mu-\alpha,\mu+\alpha)$, allowing you to account for your other sources of noise.
