Propagate errors in measured points to Simpson's numerical integral

Question

I have a set of measured/observed $y(x)$ points, each with an assigned standard deviation:

$$y: \{y_0\pm\sigma_{y_0}, y_2\pm\sigma_{y_2}, ..., y_N\pm\sigma_{y_N}\}$$

I use scipy's implementation of Simpson's rule to obtain the integral along the entire range where $y$ is sampled, $x: \{x_0, x_2, ..., x_N\}$, i.e.:

$$A = \int_{x_0}^{x_N} y(x)\,dx$$

I'd like to use the known errors in $y(x)$ to estimate a standard deviation for the integral value: $A\pm\sigma_A$.

I see in WP that the error in Simpson's rule integral can be obtained via:

$$\tfrac{1}{90} \left(\tfrac{x_N-x_0}{2}\right)^5 \left|f^{(4)}(\xi)\right|$$

where $f(x)$ is the function being integrated (which I don't have) and $\xi$ is "some number between $x_0$ and $x_N$".

I don't really understand the above expression$^{(1)}$, but as far as I can tell that's the error associated to the numerical integration, not the error given by the uncertainties in the measured points.

How can I obtain the error estimation for the integral, using the information carried by the errors in the measured points (i.e.: $\sigma_{y_i}$)?

---

$(1)$: Since I don't have $f(x)$, how am I supposed to obtain its fourth derivative? And what is $\xi$?

Answer 1

To spell out the implications of a comment by @whuber, if the $x_i$ are evenly spaced and $n$ is even, Simpson's rule estimates the value of the integral as:

$$\int_{x_0}^{x_n} y(x) \, dx \approx \tfrac{\Delta x}{3}\left(y_0 + 4y_1+2y_2+ 4y_3+2y_4+\cdots+4y_{n-1} + y_{n}\right),$$ where $\Delta x = \frac{{x_n}-{x_0}}{n}$ and $y_i=y(x_i)=y({x_0}+i\Delta x)$. If the $y_i$ are uncorrelated and have standard deviations as stated in the question, then the estimated variance of this weighted sum is:

$$ \tfrac{\Delta x^2}{9} \left(\sigma_{y_0}^2 + 16\sigma_{y_1}^2+4\sigma_{y_2}^2+ 16\sigma_{y_3}^2+4\sigma_{y_4}^2+\cdots+16\sigma_{y_{n-1}}^2 + \sigma_{n}^2\right).$$

This provides the contribution to error arising from the uncertainty in the $y_i$ values. It does not incorporate the error from the Simpson's rule approximation itself, for which you do need to know the actual functional form of $y(x)$, as noted in the question. Without that information about the form of $y(x)$, or at least some assumptions about it, I don't see how one can go farther in estimating the error.

The Simpson's rule approximation unevenly weights the variances of the $y_i$. Other approximations weight the variances more evenly. For example, with the trapezoidal rule and evenly spaced $x_i$ the corresponding approximation to the integral is:

$$\int_{x_0}^{x_n} y(x) \, dx \approx \tfrac{\Delta x}{2}\left(y_0 + 2y_1+2y_2+ 2y_3+2y_4+\cdots+2y_{n-1} + y_n\right)$$

and the variance of the approximation arising from errors in uncorrelated $y_i$ is:

$$ \tfrac{\Delta x^2}{4} \left(\sigma_{y_0}^2 + 4\sigma_{y_1}^2+4\sigma_{y_2}^2+ 4\sigma_{y_3}^2+4\sigma_{y_4}^2+\cdots+4\sigma_{y_{n-1}}^2 + \sigma_{n}^2\right).$$

The uncertainty arising from the trapezoidal rule approximation itself poses similar problems as for Simpson's rule.

If the $x_i$ are unevenly spaced, the trapezoidal approximation to the integral is:

$$\int_{x_0}^{x_n} y(x)\, dx \approx \sum_{i=1}^n \frac{y_{i-1} + y_i}{2} \Delta x_i,$$

where $\Delta x_i = x_i - x_{i-1}$. The corresponding contribution of uncertainty in $y_i$ to the variance of this estimate of the integral (with uncorrelated $y_i$) is:

$$\frac{1}{4} \left(\sigma_{y_0}^2 \Delta x_1^2+\sigma_{y_n}^2 \Delta x_n^2+\sum_{i=1}^{n-1} \sigma_{y_i}^2 (\Delta x_i + \Delta x_{i+1})^2\right),$$

so variances of $y$ values within longer $\Delta x$ intervals are more highly weighted. The paged linked above has a formula for Simpson's rule with uneven $\Delta x_i$, but the corresponding formula for the variance is left as an exercise for the reader.