Examples for error propagation

Hear we use the same config in particle_amplitude.py

config_str = """
decay:
    A:
       - [R1, B]
       - [R2, C]
       - [R3, D]
    R1: [C, D]
    R2: [B, D]
    R3: [B, C]

particle:
    $top:
       A: { mass: 1.86, J: 0, P: -1}
    $finals:
       B: { mass: 0.494, J: 0, P: -1}
       C: { mass: 0.139, J: 0, P: -1}
       D: { mass: 0.139, J: 0, P: -1}
    R1: [ R1_a, R1_b ]
    R1_a: { mass: 0.7, width: 0.05, J: 1, P: -1}
    R1_b: { mass: 0.5, width: 0.05, J: 0, P: +1}
    R2: { mass: 0.824, width: 0.05, J: 0, P: +1}
    R3: { mass: 0.824, width: 0.05, J: 0, P: +1}

"""

import matplotlib.pyplot as plt
import yaml

from tf_pwa.config_loader import ConfigLoader
from tf_pwa.histogram import Hist1D

config = ConfigLoader(yaml.full_load(config_str))
input_params = {
    "A->R1_a.BR1_a->C.D_total_0r": 6.0,
    "A->R1_b.BR1_b->C.D_total_0r": 1.0,
    "A->R2.CR2->B.D_total_0r": 2.0,
    "A->R3.DR3->B.C_total_0r": 1.0,
}
config.set_params(input_params)

data = config.generate_toy(1000)
phsp = config.generate_phsp(10000)

8.8%[▓▓▓▓>---------------------------------------------] 0.43/4.94s eff: 90.000000%
93.2%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓>---] 1.36/1.46s eff: 7.277187%
100.1%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓>] 1.80/1.80s eff: 6.216684%
100.0%[▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓] 1.80/1.80s  eff: 6.180988%

After we calculated the parameters error, we will have an error matrix config.inv_he (using the inverse hessain). It is possible to save such matrix directly by numpy.save and to load it by numpy.load.

config.get_params_error(data=[data], phsp=[phsp])

Using Model
Time for calculating errors: 0.5375051498413086
/home/docs/checkouts/readthedocs.org/user_builds/tf-pwa/checkouts/stable/docs/../tf_pwa/utils.py:265: UserWarning: matrix is not positive definited
  warnings.warn("matrix is not positive definited")
eigvalues:  [ 0.03632311  0.01773997  0.01428318 -0.00113541  0.00060804  0.00241579]
DType class corresponding to the scalar type and dtype of the same name.

Please see `numpy.dtype` for the typical way to create
dtype instances and :ref:`arrays.dtypes` for additional
information.
hesse_error: [0.05716870899070608, 0.11223857683598046, 0.10813562180020414, 0.10252326365573486, 0.07774899400317759, 0.1616158633065301]

{'A->R1_b.BR1_b->C.D_total_0r': 0.05716870899070608, 'A->R1_b.BR1_b->C.D_total_0i': 0.11223857683598046, 'A->R2.CR2->B.D_total_0r': 0.10813562180020414, 'A->R2.CR2->B.D_total_0i': 0.10252326365573486, 'A->R3.DR3->B.C_total_0r': 0.07774899400317759, 'A->R3.DR3->B.C_total_0i': 0.1616158633065301}

We can use the following method to profamance the error propagation

\[\sigma_{f} = \sqrt{\frac{\partial f}{\partial x_i} V_{ij} \frac{\partial f}{\partial x_j }}\]

by adding some calculation here. We need to use tensorflow functions instead of those of math or numpy.

import tensorflow as tf

with config.params_trans() as pt:
    a2_r = pt["A->R2.CR2->B.D_total_0r"]
    a2_i = pt["A->R2.CR2->B.D_total_0r"]
    a2_x = a2_r * tf.cos(a2_i)

And then we can calculate the error we needed as

print(a2_x.numpy(), pt.get_error(a2_x).numpy())

-0.8322936730942848 0.2416551822324546

Uncertainties of fit fractions

We can also calculate some more complex examples, such as the fit fractions of all in C+D. Even further, we can get the error of error in the meaning of error propagation.

amp = config.get_amplitude()

with config.params_trans() as pt1:
    with config.params_trans() as pt:
        int_mc = tf.reduce_sum(amp(phsp))
        with amp.temp_used_res(["R1_a", "R1_b"]):
            part_int_mc = tf.reduce_sum(amp(phsp))
        ratio = part_int_mc / int_mc
    error = pt.get_error(ratio)

print(ratio.numpy(), "+/-", error.numpy())
print(error.numpy(), "+/-", pt1.get_error(error).numpy())

0.42520552652768934 +/- 0.011657244750159976
0.011657244750159976 +/- 0.0008577773017336707

For large data size it would be some problem named OOM (out of memory). TFPWA provide vm.batch_sum_var to do sum of large samples

int_mc_v = config.vm.batch_sum_var(amp, phsp, batch=5000)

with amp.temp_used_res(["R1_a", "R1_b"]):
    part_int_mc_v = config.vm.batch_sum_var(amp, phsp, batch=5000)

It will store the pre-calculated gradients as

print(int_mc_v.grad, part_int_mc_v.grad)

tf.Tensor(
[1751674.66853607 -117780.40351603 1359484.69263483   79531.56972059
  783486.71769585   10854.12939637], shape=(6,), dtype=float64) tf.Tensor(
[1848665.96791494   -4822.61854968       0.               0.
       0.               0.        ], shape=(6,), dtype=float64)

Then, we can use it as a function to do error propagation:

with config.params_trans() as pt:
    ratio = part_int_mc_v() / int_mc_v()
error = pt.get_error(ratio)

print(ratio.numpy(), "+/-", error.numpy())

0.42520552652768934 +/- 0.011657244750159973

Besides the error propagation, there would be some additional uncertainties. For example, the uncertainty from the integration sample size is often defined as the sum of square as

with amp.temp_used_res(["R1_a", "R1_b"]):
    int_square = tf.reduce_sum((amp(phsp) / int_mc) ** 2)

print(ratio.numpy(), "+/-", error.numpy(), "+/-", tf.sqrt(int_square).numpy())

0.42520552652768934 +/- 0.011657244750159973 +/- 0.010392777477148951