Uncertainty Propagation according to GUM

2

The difference of two Poisson RVs is distributed by the Skellam distribution. We have an implementation in pymc-experimental which you could directly use to model your data. Something like:

import pymc as pm
import pymc_experimental as pmx

 with pm.Model() as m:
    mu1 = pm.HalfNormal('mu1', 1)
    mu2 = pm.HalfNormal('mu2', 1)
    y_hat = pmx.distributions.Skellam('y_hat', mu1, mu2, observed=data)