# Error propagation with an integral

My question here is about how to determine the error of an integral given individual uncertainties in two parameters defining the function being integrated.

I used a curve fitting function to approximate a peak in my dataset as a given function that is essentially a Gaussian. This has given me uncertainties in both the height and width of the function. I have then integrated the fit function, and I would like to find an uncertainty associated with the value of this integral. How would I go about doing this given the uncertainties for height and width of my peak?

• I'm likely to use the above mentioned method as well, and I'd REALLY like to see a good answer to this question! Jul 31, 2015 at 3:45
• There is no obvious answer to this, for, say, $\int^{x_{1}}_{x_{0}} (f \pm g)dx$ for general error function $g$. if the function is monotonic the associated error is easy and will just have maximum error $G(x_{1})$ with $G = \int g$ and similar for the minimum. If the function is complicated you will have to determine maxima and minima and such. Jul 31, 2015 at 3:46
• @DavidWhite You know you can mark a question as a favorite, don't you, to come back to? Jul 31, 2015 at 6:29

Take a Gaussian integral $$I(a,b) = \int^\infty_{-\infty} a \,e^{-(\frac{x}{b})^2} dx \,= ab\sqrt{\pi}$$ It is a function of parameters $a,b$ that have some error $\sigma_a,\sigma_b$. In principle the errors might have some correlation $\sigma_{ab}$ although this is considered to be zero in curve fitting as far as I know. The usual expression for the error $\sigma_I$ of a function of parameters with error is $$\sigma_I^2 = \left(\frac{\partial I}{\partial a}\right)^2 \sigma_a^2 + \left(\frac{\partial I}{\partial b}\right)^2 \sigma_b^2 + 2 \frac{\partial I}{\partial a}\frac{\partial I}{\partial b}\sigma_{ab}$$ In the case of the Gaussian integral with uncorrelated error in the parameters $$\sigma_I = \sqrt{\pi(b^2\sigma_a^2 + a^2\sigma_b^2)}$$

• Just for curiosity's sake, where did the above formula come from? Jul 31, 2015 at 4:19
• @David White, Expand $I(a,b)$ as a taylor series keeping only the first order terms. Subtract the mean value $I(\bar{a},\bar{b})$ then square and take the expectation value (the definition of the variance). Jul 31, 2015 at 4:29

The problem with Octonion's answer is that you don't in general know the full covariance matrix for the errors: you can readily estimate $\sigma_a$ and $\sigma_b$ but $\sigma_{a\,b}$ is a little harder. You can handle this by overbounding: by the Cauchy-Schwarz inequality, we can say $\sigma_{a\,b} \leq \max(\sigma_a,\,\sigma_b)^2$ and then you replace Octonion's formula by:

$$\sigma_I\leq\sqrt{\pi(b^2\,\sigma_a^2 + a^2\, \sigma_b^2 + 2 |a|\,|b|\,\max(\sigma_a,\,\sigma_b)^2}$$

but this may be too pessimistic for you, especially if there is a large negative correlation between the two estimated parameters.

An alternative approach is linear regression: you may be able to represent your function with a superposition of, say, $M$ Tschebyschev polynomials or other orthogonal set like a Fourier series and fitting this superposition to $N$ datapoints. Then, let $\Phi$ be the $N\times M$ matrix such that $\Phi_{j\,k}$ is the value of the $k^{th}$ basis function at the $j^{th}$ measurement point then:

$$\Phi\,\alpha = \hat{Y}$$

where $\alpha$ is the $M\times 1$ column vector of superposition weights (regression co-efficients) for the $M$ basis functions and $\hat{Y}$ the $N\times 1$ column vector of estimated values for the experimental measurements. Then if $Y$ is the $N\times 1$ column vector of measurements, we seek to minimize:

$$\epsilon^2 = \|\Phi\,\alpha-Y\|^2 = (\Phi\,\alpha-Y)^T\,(\Phi\,\alpha-Y)$$

and given

$$\mathrm{d}_\alpha \epsilon = \alpha^T\,\Phi^T\,\Phi - Y^T\,\Phi$$

and we seek an extremal value of $\epsilon^2$, so the value of $\alpha$ that yields this fulfils:

$$\Phi^T\,\Phi\,\alpha = \Phi^T\,Y$$

and you get your superposition weights ("regression co-effcients") by:

$$\alpha = (\Phi^T\,\Phi)^{-1}\,\Phi^T\,Y$$

You need to choose enough measurement points and the right basis functions so that $\Phi^T\,\Phi$ is well conditioned. If this matrix is not well conditioned, it means either that you don't have enough data points OR that your basis functions are only weakly determined by your measurements OR that your measurement points are positioned badly. Example: any constant times a basis function $\phi_j(x)$ which is nought at all the measurement points adds nothing to the regression error: conversely, nothing can be inferred about such a function's weight in the superposition from the particular measurement points in question.

The way to look at the error behavior of this equation is with singular value decomposition of both sides which writes:

$$(\Phi^T\,\Phi)^{-1}\,\Phi^T = U\,\Sigma\,V$$

where $U,\,V$ are $M\times M$ orthonormal (real unitary) matrices, and $\Sigma$ a diagonal matrix of nonzero, positive singular values. Note that, beginning with the assumption that the measurement errors are uncorrelated, the matrices $U,\,V$ are unitary, which means that they map identical uncorrelated random variables to indentical uncorrelated random variables. So now we have $U^T\,\alpha$ is a set of uncorrelated RVs with variances $\sigma_1^2\,s^2,\,\sigma_2^2\,s^2,\,\cdots$ where the $\sigma_j^2$ are the $M$ singular values in $\Sigma$ and $s^2$ is the variance of your measurement errors.

So now let $\phi_j(x)$ be the basis functions. Then the integral you need is $I^T\,\alpha$, where $I(x)$ is the $M\times 1$ column vector of integrals $I_j = \int_{x_1}^{x_2}\phi_j(u)\mathrm{d}u$. So now work out $I^T\,U$: this will be a row vector $(v_1,\,v_2,\,\cdots)$. When it foremultiplies $U^T\,\alpha$ it gives you $I^T\,\alpha$, which is the integral you need. So the variance of your integral is now:

$$\sigma_I^2 = (v_1^2\,\sigma_1^2 + v_2^2\,\sigma_2^2 + \cdots)\,s^2$$

where the $\sigma_j^2$ are the $j$ singular values, $s^2$ is the variance of the error in each measurement point and the $v_j$ are the members of the row vector $I^T\,U$.

If you are happy with a numerical recipe then by far the most genuine, easy and error-less method is a pseudo-data approach.

Theoretical approach is very often just an approximation (linearization in formulas) and cannot take into account full distributions. Even full correlation matrix contains only a very concentrated information from situation which is characterized by smooth functions (smooth probability distributions).

Recipe is simple: generate you parameters randomly with appropriate means and standard deviations. For each pair (if you have two of them) compute integral $$I_j$$. After $$N$$ evaluations ($$N$$ times you generate both parameters and compute integral) you can compute the mean and the error of your integral:

$$I_{mean} = \frac{1}{N} \sum_{j=1}^{j=N} I_j$$ $$I_{error} = \sqrt{ \frac{1}{N} \sum_{j=1}^{j=N} (I_{mean} - I_j)^2 }$$

Of course you should try to make $$N$$ big. With $$N$$ large enough, this approach provides you with the full $$I$$ distribution as function of distributions of you two parameters. This approach, of course, takes any possible correlation of your input parameters into account correctly and automatically (if you have correlated input parameters you should generate them in a correlated way).

Simply said: because this approach stems from the very definition of the error propagation itself, its only limitation is statistics (or computer time).