How does uncertainty of dataset propagates through numerical integration? 
In the following text, uncertainty refers to standard deviation.

I have 500 time series which I use in a few equations and get averages and uncertainties. Through calculations for each time series, I end up with a data set with 500 values of the position z, the resistance $R(z)$, the absolute and the relative uncertainty:
$$  z \,,  R( z )  \,, \pm \Delta R(z) \,, \pm \frac{\Delta R(z)}{ R(z) }$$
What I want is to calculate the permeation $P$ through the following formula:
$$ P=\frac{1}{\int_{z_{min}}^{z_{max}}R(z)dz} $$
therefore, I need to numerically integrate these 500 values. This is easy with a trapezoidal but the question is what happens to the uncertainty? How does it propagate through the integral? 
Essentially I want to report the average $ \left( P \pm \Delta P \right) \, [units]$
 A: 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 ab 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)

