# 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?

• Just clarifying: do you mean that there is no error other than the uniform error, and you know the bounds on it? There's no possibility of the true curve going outside the interval? – jbowman May 20 '12 at 12:58
• I could also have another source of noise which means that I could have the true curve appearing outside the interval. – alang May 20 '12 at 14:14

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.

• Thanks for the answer. Your Gaussian process formulation seems interesting, but as you say, it does not seem so easy to compute. The quadratic programming approach does not seem to extend easily in the presence of other sources of noise, unless I'm mistaken. – alang May 22 '12 at 7:49
• Yes, the QP solution is only for the uniform case, but this is also true for the the described Gaussian Process approach. When its dispersion is large, a normal truncated to $(y_{i0},\,y_{i1})$ becomes nearly uniform, and this is what happens in the (constrained) Smoothing Spline when $p$ is small. A possible formulation for the two sources of noise in your problem is using a functional prior for $\eta(x)$, as well as a distribution for the observation(s) at $x$ conditional on $\eta(x)$, typically $\eta(x)+$ noise. – Yves May 22 '12 at 8:39

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.

• Kernel smoothing is an option, but it's not clear to me how to build a kernel which will guarantee the property I want without being too restrictive. I think I would also run into the same problem of favouring the mean of the interval. – alang May 20 '12 at 14:10

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.

• This seems like it would work well, but is there any theoretical justification for using a generalized normal distribution? – alang May 21 '12 at 8:31
• It's a simple, well-behaved function that approximates your error model with nonzero density throughout the real line, so it will not fail to locate points outside the uniformly-distributed range. A concern is that it may not give these aberrant points their due weight, but we need more information from you to fix that. – Emre May 21 '12 at 16:59
• Suppose I know that the extra noise is normally distributed and I have an estimator for its variance, does that determine the parameter $\beta$? – alang May 22 '12 at 5:47
• Yes, it does. Assuming your noise terms are independent, their variances are additive. Furthermore you have analytic expressions for the variances of all the random variables. That of the generalized normal distribution depends on gamma through the gamma function, so you might have to use a numerical search. Alternatively, if you know your extra noise is normally distributed, you could simply add another noise term. – Emre May 22 '12 at 6:00
• I'm not sure I understand your last comment. In the case that my noise is normally distributed, you mean that I can add a noise term in addition to the generalized normal distribution? It's still not clear to me how to choose the parameters $\alpha,\beta$ for the generalized normal distribution. – alang May 22 '12 at 6:12