Klaus Reygers, 2020
from sympy import *
from IPython.display import display, Latex
def gaussian_error_propagation_corr(f, x, V):
"""
f: function f = f(x[0], x[1], ...)
x: list of variables
V: covariance matrix (python 2d list)
"""
sum = S(0) # empty sympy expression
for i in range(len(x)):
for j in range(len(x)):
sum += diff(f, x[i]) * diff(f, x[j]) * V[i][j]
return sqrt(simplify(sum))
Show usage for a simple example: Volume of a cylinder with radius $r$ and height $h$:
r, h, sigma_r, sigma_h = symbols('r, h, sigma_r, sigma_h', positive=True)
rho = Symbol("rho", real=True) # correlation coefficient
V = pi * r**2 * h # volume of a cylinder
vars = [r, h]
cov_matrix = [[sigma_r**2, rho * sigma_r * sigma_h],
[rho * sigma_r * sigma_h, sigma_h**2]]
Matrix(cov_matrix)
sigma_V = gaussian_error_propagation_corr(V, vars, cov_matrix)
display(Latex(f"$V = {latex(V)}, \, \sigma_V = {latex(sigma_V)}$"))
Plug in some numbers and print the calculated volume with its uncertaity:
r_meas = 3 # cm
sigma_r_meas = 0.1 # cm
h_meas = 5 # cm
sigma_h_meas = 0.1 # cm
central_value = V.subs([(r,r_meas), (h, h_meas)]).evalf()
sigma = sigma_V.subs([(r, r_meas), (sigma_r, sigma_r_meas), (h, h_meas), (sigma_h, sigma_h_meas)]).evalf()
for rho_value in [-1, 0, 1]:
sigma = sigma_V.subs([(r, r_meas), (sigma_r, sigma_r_meas), (h, h_meas), (sigma_h, sigma_h_meas), (rho, rho_value)]).evalf()
display(Latex(f"$$ \\rho = {rho_value}: V = ({central_value:0.1f} \pm {sigma:.1f}) \, \mathrm{{cm}}^3$$"))