Numerical integration with measurement errors: what scheme to use?

Question

I have an unknown 'smooth' function that is sampled at some $n$ equally spaced locations $x_i$ between $0$ and $1$, with some measurement values $f(x_i)=y_i$, but each has some Gaussian uncertainty ($(\sigma_y)_i$). I am interested in the integral $I=\int f\; \mathrm{d}x$.

One way to get it is to use a Riemann sum, which would yield $I=\frac{1}{N}\sum_i x_i$. (Or if the errors are different, the sum should be weighted.) But Riemann integration just assumes a piecewise function which is not a good approximation if there is some curvature in the function.

Another way would be to use a higher order method, for example Simpson's rule. These are higher order, but have worse behaviour with large uncertainties: if one measurement is off, the quadratic interpolation would be much more wrong than a linear interpolation, and hence the integral of the quadratic interpolation would be worse than the integral of the linear interpolation.

So there seems to be a tradeoff: if there are high uncertainties, use lower-order methods, and if the uncertainties are low, use higher order methods. But how to choose? What would be the best way of estimating $I$ (and its uncertainty)?

One way would be to parametrize the 'smooth function', estimate the parameters and calculate the integral of the parametrization. But then, if there is not an obvious model, what parametrization should I use? (And what priors on the parameters?)

Answer 1

I think what you need is the smoothing spline. You'd be minimizing the following functional by looking for a function $\hat f(x)$: $$ \sum_{i=1}^n \{y_i - \hat f(x_i)\}^2 + \lambda \int \hat f''(x)^2 \,dx. $$ Once you got the function, integrating it is yrivial for cubic or other polynomial splines.

There a similar technique called Kernel density estimation. Here you approximate your function with a chosen kernel:
$$ \widehat{f}_h(x) = \frac{1}{n}\sum_{i=1}^n K_h (x - x_i) = \frac{1}{nh} \sum_{i=1}^n K\Big(\frac{x-x_i}{h}\Big), $$ Again, once you have a Kernel, integrating it is trivial.

Neither of these techniques use the distribution of errors directly, i.e. you don't need them to be Gaussian.