Uncertainty calculation through integration and correct analysis methodology

Question

In the following text, uncertainty refers to standard deviation. My methodology and notation was based mainly on the book "A Student's Guide to Data and Error Analysis" by H.J.C. Berendsen."

Summary:

I have a timeseries which I use in a few equations and get averages and uncertainties. I do that up to a point where I have to integrate the averages. So the question is, what happens to uncertainty when you integrate the averages? My problem might arise from a mistake in methodology so I would appreciate any feedback on the followed process as well.

Details:

I am trying to calculate the uncertainty of a number that is computed as part of an analysis code. While I can compute the average number, I encounter some issues with the uncertainty for which I would appreciate any advice from more experienced analysts.

My methodology was:

1. From 30 simulations, I have obtained 30 timeseries files that each contain 100000 lines of a z-axis position-fluctuation of a molecule. Therefore I have 30 files files that each contains two vectors: t and z

2. For each position/file, I use the z timeseries to calculate  the variance of the timeseries $var(z)$ and the variance of the mean of the timeseries $var(\langle z \rangle)$. The latter is done with block averaging (b.a.). Both these methods have a textbook defined, relative standard deviation: $$ var(z) \, [units] \pm \sqrt{\frac{2}{n-1}} $$ $$ var(\langle z \rangle) \, [units] \pm \sigma_{<χ>}^2 $$ where $ \sigma_{<χ>}^2 $ is the squared standard error ($ \sigma_{<χ>}=\frac{\sigma}{\sqrt{n}}$) computed from b.a.  The above are used as input in a diffusion equation: $$ D(z)=\frac{var(z)}{\big[ \frac{n \cdot var(\langle z \rangle) }{var(z)} -1 \big] \frac{dt}{2}} $$ where $n$ and $dt$ are constants. For this step, I also used two failed alternative methods that I explain in the end.

3. Then I used the above to calculate the average diffusion coefficient D for each of the 30 files. To find the uncertainty of each average, I combined the errors based on error propagation algebra where I considered (as a first approximation) the covariance $\sigma_{AB}$ as 0. In the end, I obtained a new "summary" file with 30 rows and 4 columns, containing the average position, average Diffusion, absolute uncertainty and relative uncertainty: $$ z \,, D(z) \,, \pm \Delta D(z) \,, \pm \frac{\Delta D(z)}{ D(z) }$$

4. From the same simulation, I also obtain a similar "summary" file that contains 30 rows and 3 columns with the average position, Gibbs free-energy and its absolute uncertainty (standard deviation): $$ z \,, \Delta G(z) \,, \pm \Delta \Delta G(z) $$

5. I then combined both summary files to compute a resistance profile with the equation: $$ R(z) = \frac{\exp \big[ {\frac{\Delta G(z)}{k_BT}}\big] }{D(z)} $$ where $k_BT$ was a constant. To combine $\Delta G(z)$ and $D(z)$ I first I interpolated them in a common, dense z-grid. Therefore, I ended up with a file with 500 values of resistance together with the absolute and relative uncertainty: $$ z \,, R( z ) \,, \pm \Delta R(z) \,, \pm \frac{\Delta R(z)}{ R(z) }$$

6. To get the final result which is permeation, I integrated these 500 values, according to the equation: $$ P=\frac{1}{\int_{z_{min}}^{z_{max}}R(z)dz} $$ and I obtained the average P by performing a trapezoidal approximation over the 500 $R(z)$ values.

My question is, what happens to the uncertainty? How does it propagate through the integral?

Also, do you have a more appropriate alternative to my methodology? I am trying to find some guidelines on how to do a "correct" data analysis but the available information is either too theoretical or too vague. For example, would it make any difference in this case to perform fitting on the $D(z)$, $\Delta G(z)$ and $R(z)$ profiles, before integration? If yes, how would I calculate the uncertainty for each and for the $P$?

Thanks a lot for any help or advice! :D

In step 2, I tried 2 more things:

1. instead of calculating the uncertainties of $var(z)$ and $var(\langle z \rangle)$ from the textbook, I used bootstrapping to create 100 new timeseries of z, for each position/file. And then calculated 100 different diffusions which I used to calculate the average and s.d. The method was returning the same statistics for the timeseries (which was expected) but the final Diffusion was extremely oscillatory (was even negative), therefore I dropped the method.

2. I tried to use block averaging for the calculation of D, but Python was refusing to do recursive block averaging (I might had a bug in my code as well). Anyway, I quickly abandoned the idea as well.

Answer 1

I found an answer to my question which in my opinion is correct (or at least has a point). Nevertheless it would be nice to have a validation from someone that might know more on the topic.

Since the average is calculated with a trapezoidal (which in reality is a summation) then the standard deviation can also be calculated by uncertainty propagation algebra for addition.

Based on Wikipedia for trapezoidal rule: $$ \int_c^d f(x) dx \approx (d-c) \left[ \frac{f(c)+f(d)}{2} \right] = \left( \frac{d-c}{2} \right) f(c) + \left( \frac{d-c}{2} \right) f(d)$$

which will follow the 2nd uncertainty propagation rule from Wikipedia: $$ f=aA+bB$$

with:

• $ a = b = \frac{d-c}{2}$
• $ A = f(c) $
• $ B = f(d) $

In terms of Python3 code that translates to:

import numpy as np

# FUNCTIONS
def uncAddition(A=1.,B=1.,sigmaA=0.,sigmaB=0.,a=1.,b=1.,sigmaAB=0.):
    f =  (a * A) + (b * B)
    sigmaf = np.sqrt( (a * sigmaA)**2 \
                    + (b * sigmaB)**2 \
                    + (2 * a * b * sigmaAB) )
    return f, sigmaf


# MAIN
# Data: is an nx2 array that contains n number of f values 
#           and their absolute standard deviation
# x:    is the independent variable of the f values
# Obviously x and f should have the same number of rows

integralAve=0.
integralStD=0.
for i in range(1,len(data[:,0])):
    trapzPartAve,trapzPartStd = uncAddition(A=f[i-1,0],
                                            sigmaA=f[i-1,1],
                                            a=(x[i]-x[i-1])/2.0,
                                            B=f[i,0],
                                            sigmaB=f[i,1],
                                            b=(x[i]-x[i-1])/2.0) 
    integralAve += trapzPartAve
    integralStd += trapzPartStd

print(integralAve,integralStd)