0
$\begingroup$

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

$\endgroup$
3
  • $\begingroup$ (1) Are the observations independent? (2) What do you mean by "chi"? There is no such symbol anywhere in your question. (3) Note that the formula for the error assumes "the points at which the integrand is evaluated are distributed symmetrically in the interval... ." Much of what you ask about Simpson's Rule is addressed in the subsequent discussion in the Wikipedia article. $\endgroup$
    – whuber
    May 26 '16 at 19:30
  • $\begingroup$ (1) yes they are, (2) I meant $\xi$, fixing it now, (3) not sure what you mean by "addressed in the subsequent discussion in the Wikipedia article". I can see nothing there about how to use the errors in the measured points to asses an error for the integral. $\endgroup$
    – Gabriel
    May 26 '16 at 19:43
  • $\begingroup$ The article addresses numerical errors arising from differences between $f$ and its approximating fourth-degree polynomial. It contains answers to what $\xi$ ("xi") means and the meaning of the fourth derivative. It is somewhat sketchy, but any good Calculus textbook will contain all the missing details and explanation. Apart from those issues, Simpson's Rule estimates the integral with a linear combination of the $y_i$, so your question comes down to computing the variance of a linear combination of independent random variables. $\endgroup$
    – whuber
    May 26 '16 at 20:01
2
$\begingroup$

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.

$\endgroup$
9
  • $\begingroup$ I have trouble understanding the final equation. It looks like this was computed for every $i$ independently. But, e.g., $y_i^2$ appears to two summands, implying 100% correlation. If I first combine the two contributions to each $y_i$ and then compute the variance I seem to arrive at a different result? $\endgroup$
    – Simon
    Dec 13 '21 at 14:38
  • $\begingroup$ Sorry, I meant to write $y_i$ (no square). Let us say we compute the contribution by $\sigma_j^2$. Since $y_j$ appears in the summands for $i=j$ and $i=j+1$, I get a term $(x_{j+1}-x_{j-1})^2 \sigma_j^2/4$. $\endgroup$
    – Simon
    Dec 13 '21 at 15:22
  • $\begingroup$ My points is: Yes, we assume $y_i$ and $y_j$ to be uncorrelated, but $y_i$ is correlated with itself, and every $y_i$ appears in two summands in the composite trapezoid rule, so we cannot compute the error for every summand independently and sum the resulting errors. $\endgroup$
    – Simon
    Dec 13 '21 at 15:24
  • 1
    $\begingroup$ @Simon you are correct. I think that I have now corrected the error. Thanks for finding it. $\endgroup$
    – EdM
    Dec 13 '21 at 23:03
  • 1
    $\begingroup$ You are right, I found my error. $\endgroup$
    – Simon
    Dec 14 '21 at 19:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.