# How to integrate numerically a function with error bars?

Typically, the function that one wants to integrate numerically, $f$, is given, i.e. its values for various points $\{x_i\}$ are known precisely. The resulting error is due to the fact that we chose a finite number of points. Let's call this error $numerical$.

I have a different situation. I know the values of $f$, $\{f(x_i)\}_{i=1}^N$, but with corresponding variances, $\{var(f(x_i))\}$. If I use, for example, the trapezoidal rule: $$F(b)-F(a)=\int_a^b f(x) dx \approx \frac{\Delta x}{2} [f(x_1) + 2 f(x_2) + \ldots 2f(x_{N-1}) + f(x_N) ]$$ (assuming that $x_i$ are equally separated by $\Delta x$), I calculate the variance of this estimate as: $$\left(\frac{\Delta x}{2}\right)^2 \left[ var(f(x_1)) + 4\cdot var(f(x_2)) + \ldots + 4\cdot var(f(x_{N-1})) + var(f(x_N)) \right]$$ but what about the $numerical$ error?

How can I combine these two errors? What is the "final" error of my estimate?

EDIT:

I've calculated the variance the same way I would do it for a sum of $N$ independent random variables, $X_i$, with coefficients $a_i$. That is, if: $$S :=\sum_{i=1}^{N}a_i X_i,\qquad \forall_{i\neq j}Cov(X_i,X_j)=0$$ then $$Var(S)=\sum_{i=1}^{N} a_i^2 Var(X_i).$$ If this error analysis is not correct, please let me know.

Also, I noticed a problem with the variance of partial sums, $S_k:=\sum_{i=1}^k a_i X_i$. It seems that the error propagates, and therefore $S_N$ has the largest variance, while $S_2$ -- the lowest. But if I were to reverse the order in the sum (which would correspond to estimating a reversed integral, $\int_b^a f dx$), the error would again propagate, but in the opposite direction. I would get two different estimates of "partial integrals", depending on the direction of integration.

Is there a way of combining information from numerical integration in both ways?

Assuming the integration is deterministic (that is, ignoring effects like bit flips and hardware issues that are generally negligible at scale), you have two types of uncertainties:

• the uncertainty in calculating the integral approximately with numerical methods instead of exactly; this error is epistemic, and with enough analysis, could be calculated or estimated

• the uncertainty in each of the function values sampled; unless there's additional relevant information, this uncertainty appears to be aleatoric

Since you have variances, one approach to quantifying the aleatoric uncertainty that is compatible with the information you've provided is polynomial chaos. You could integrate the mean and variance separately to obtain the mean and variance of the integral. This approach is different than the one you propose above; how you arrived at the formula above for the variance is unclear to me.

Another approach would use interval arithmetic, but this approach is better at capturing worst-case upper and lower limits, and it's not clear that such information would be valuable to you.

As for the remaining epistemic uncertainty, you could ignore it, if it's sufficiently small relative to the aleatoric uncertainty in the integral (and in your data points). If you feel you do need to quantify it, you'd need to estimate or bound the error somehow (estimates via embedded methods, arbitrary precision integration, or bounds from numerical analysis textbooks/research). As a first approximation, you could model the epistemic error in the integral of the mean as a uniform random variable at each point in space, calculate the resulting variance, and combine that variance information with the integral of the variance calculated from the aleatoric uncertainties, using a polynomial chaos approach. The error in integrating the variances could then be incorporated into calculations of higher-order moments of the random variable describing the integral you wish to calculate. Although this process could be continued with progressively higher-order moments, at some point, the moment expansion of this variable will need to be truncated in calculations.

• I'll take a look at polynomial chaos. In the mean time, I've corrected an error in the variance formula (replaced "4" by "2" in the fraction). Dec 4 '14 at 8:10
• Depending on the application and the size, couldn't it end up being simpler to just perform some Monte Carlo sampling? i.e. just perform the numerical integration many times, each time using values of $f$ selected from the distribution? I haven't played with polynomial chaos, but it seems like it could be overkill in many cases.
– AJK
Dec 6 '14 at 3:42
• @AJK: It could be. Monte Carlo tends to have error proportional to the square root of the number of samples, so the number of samples for each knot (sample point) needed to sample acceptable error could be large (10,000 to 1 million). This estimate assumes a naive approach; you might be able to use variance reduction methods to use fewer samples and achieve a more accurate result, but that's outside my area of expertise. Dec 8 '14 at 9:19
• @AJK: I haven't thought about Monte Carlo, but I did consider bootstrapping. Alas, I didn't find any quick-fix solution. For now, I'm going to stick to the simple error-analysis that I've mentioned above, and later expand the method either using polynomial chaos approach, or bootstrapping. Dec 8 '14 at 9:58