Propagation of error for an integral using Leibnitz' rule

Question

I'm computing integrals $u_i$, where $$u_i = \int_a^b s^i G(s) ds$$ and $a$ and $b$ are constants.

Now I want to estimate the error in $u_i$; I have estimates for the error $\sigma_G$ of $G$, and there is no error in $s$.

Using Leibniz' rule for differentiation under the integral $${du_i \over dG} = \int_a^b {\partial(s^iG)\over \partial G}ds$$ and since $\displaystyle {\partial(s^iG)\over \partial G}=s^i$, $${du_i \over dG} = \int_a^b s^i ds = \ln(b)-\ln(a){\rm \ \ if\ }i=-1{\rm,\ else\ }{1\over i+1}(b^{i+1}-a^{i+1})$$ I want to apply the propagation of error approximation $$\sigma_{u_i}^2\approx \sigma^2_G\left({\partial u_i \over \partial G}\right)^2 $$ so I used $\displaystyle {du_i \over dG}$ for $\displaystyle{\partial u_i \over \partial G}$. I think there is something wrong with this derivation, since numerical tests do not agree with the result calculated from the above. Can anyone see why this derivation is wrong, or inapplicable?

Edit: The form of Leibnitz' rule I am using is $${d\over dx}\left(\int_{f(x)}^{g(x)} h(x,t)\;dt\right) = g^\prime(x)h\bigl(x,g(x)\bigr)-f^\prime(x)h(x,f(x))+\int_{f(x)}^{g(x)} {\partial h\over \partial x}(x,t)\;dt$$ Since $f(x)$ and $g(x)$ are the constants $a$ and $b$, all but the final term on the RHS vanish. I don't see any ambiguity in the meaning of the partial derivative in this term.

Answer 1

Here I will consider evaluating the integral by discretizing at prescribed intervals according to the information provided in an edit, considering the values of $G$ as Gaussian random variables and computing the variance of the resulting sum as a function of the number of points in the interval $N$.

Assuming that $G(s)$ is distributed as noted above, we note that with the discretization presented in the comments the integral is be approximated by the quantity

$$ I_N(a,b)=\frac{1}{2}\sum_{k=0}^{N-1}\left[(x_k)^n G(x_k)+(x_{k+1})^n G(x_{k+1})\right](x_{k+1}-x_k)$$

with $x_k=a(b/a)^{k/N}$. This can be rewritten explicitly as

$$ I_N(a,b)=\frac{1}{2}\left(a^{n+1} G(a)\left(\left(b/a \right)^{1/N}-1 \right)+b^{n+1}\left(-\left(b/a\right)^{-1/N}+1\right) G(b)+\sum_{k=1}^{N-1}G(x_k)(a(b/a)^{k/N})^{n+1}((b/a)^{1/N}-(b/a)^{-1/N})\right)$$

We now use the result that the variable $U=\sum_k r_k G_k$ is Gaussian with variance $\sum_k(\sigma_{G_k}^2r_k^2)$ to compute the variance

$$\frac{\sigma^2_{I_N}}{\sigma_G^2}=\frac{a^{2n+2}}{4}\left(\left(\frac{b}{a}\right)^{\frac{1}{N}}-1\right)^2+\frac{b^{2n+2}}{4}\left(-\left(\frac{b}{a}\right)^{-\frac{1}{N}}+1\right)^2+\frac{1}{4}\left(\left(\frac{b}{a}\right)^{1/N}-\left(\frac{b}{a}\right)^{-1/N}\right)^2\left[\frac{b^{2n+2}-a^{2n+2}}{(b/a)^{2n+2/N}-1}-a^{2n+2}\right]$$

For large $N$, we find that the variance tends to zero, with leading contribution

$$\frac{\sigma^2_{I_N}}{\sigma_G^2}=\frac{\ln(b/a)}{N}\left(\frac{b^{2n+2}-a^{2n+2}}{2n+2}\right)+\mathcal{O}(1/N^2)$$

The last two expressions are only valid for $n\neq - 1$. In order to obtain a valid formula for $n=-1$, we take the limit $n\to - 1$ of the above and we find

$$\frac{\sigma^2_{I_N}}{\sigma_G^2}=\frac{1}{4}\left(\left(\frac{b}{a}\right)^{\frac{1}{N}}-1\right)^2+\frac{1}{4}\left(-\left(\frac{b}{a}\right)^{-\frac{1}{N}}+1\right)^2+\frac{N-1}{4}\left(\left(\frac{b}{a}\right)^{1/N}-\left(\frac{b}{a}\right)^{-1/N}\right)^2$$

and in the large $N$ limit

$$\frac{\sigma^2_{I_N}}{\sigma_G^2}=\frac{\ln^2 b/a}{N}+\mathcal{O}(1/N^2)$$

We observe that:

1. The variance tending to zero as $1/N$ is justified by a slight generalization of the Central Limit Theorem. The effect should persist even for non-Gaussian variables.

2. It is impossible to write the variance simply as an integral- note that the sum contains a term that is quadratic in the spacing $\sim\sum(x_{k+1}-x_k)^2$, which is not directly amenable to a Riemann integral treatment. In fact, this object has completely different properties, in that it's explicit value depends on the discretization scheme chosen. It's an interesting exercise to partition the interval linearly instead of logarithmically to see that it produces a different result.

3. Both results for finite and large $N$ agree very well with the numerical results cited below in the comments.

Answer 2

Here is my derivation for the propagation-of-error method. Eqn (3) is the trapezoidal rule. There is no restriction on the spacing of $s$. $$\eqalignno{ \mu_i &\approx \sum_{j=1}^{n-1} {\textstyle{1\over 2}}( s^i_j G_j + s^i_{j+1} G_{j+1}) ( s_{j+1}- s_j) &(3)\cr &= {\textstyle{1\over 2}} s^i_1 G_1 ( s_2- s_1) + {\textstyle{1\over 2}} \left(\sum_{j=2}^{n-1} s^i_j G_j ( s_{j+1}- s_{j-1})\right) + {\textstyle{1\over 2}} s^i_n G_n ( s_n- s_{n-1})&(17a)\cr &= \sum_{j=1}^n a_j G_j&(17b)\cr }$$ where the coefficients $a_j$ are independent of $G_j$, and defined by $$a_j = \cases{ {\textstyle{1\over 2}} s^i_1 ( s_2- s_1) &if $j=1$\cr \cr {\textstyle{1\over 2}} s^i_n ( s_n- s_{n-1}) &if $j=n$\cr \cr {\textstyle{1\over 2}} s^i_j ( s_{j+1}- s_{j-1}) &otherwise\cr }$$ So, since ${\displaystyle\partial(a_jG_j)\over \displaystyle\partial G_j} = a_j$, $$\eqalignno{ \sigma^2_{\mu_i} &= \sum_1^n a_j^2\sigma_{G_j}^2 & (19) }$$