three ways to propagate yields

yield_prop.py

import json
import pyhf
from pyhf.contrib.utils import download
import cabinetry

download("https://www.hepdata.net/record/resource/1935437?view=true", "bottom-squarks")
ws = pyhf.Workspace(json.load(open("bottom-squarks/RegionC/BkgOnly.json")))
patchset = pyhf.PatchSet(json.load(open("bottom-squarks/RegionC/patchset.json")))
ws = patchset.apply(ws, "sbottom_600_280_150")
cabinetry.workspace.save(ws, "bottom-squarks.json")


import jax
import jax.numpy as jnp
import json
import pathlib
import jacobi
import numpy as np

# get statistcal model + data
fname = pathlib.Path("bottom-squarks.json")
spec = json.loads(fname.read_text())
ws = pyhf.Workspace(spec)
model = ws.model()
data = ws.data(model)
# fit with pyhf
pyhf.set_backend("numpy", "minuit")
result, result_obj = pyhf.infer.mle.fit(data, model, return_result_obj=True)
# error propagation
y, ycov = jacobi.propagate(
 lambda p: model.expected_data(p, include_auxdata=False),
 result_obj.minuit.values,
 result_obj.minuit.covariance,
)
print(f"via error propagation:\nyield: {y}\nunc: {np.diag(ycov)** 0.5}\n")
# bootstrap sampling
rng = np.random.default_rng(1)
par_b = rng.multivariate_normal(
 result_obj.minuit.values, result_obj.minuit.covariance, size=50000
)
y_b = [model.expected_data(p, include_auxdata=False) for p in par_b]
yerr_boot = np.std(y_b, axis=0)
print(f"via bootstrapping:\nyield: {np.mean(y_b, axis=0)}\nunc: {yerr_boot}")


pyhf.set_backend("jax")

def fisher_cov(model, pars, data):
    return jnp.linalg.inv(-jax.hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data))

F_a = jax.jacrev(lambda pars: model.expected_data(pars, include_auxdata=False))(result)
cov = fisher_cov(model, result, data)
error_squared = F_a @ cov @ F_a.T
error = jnp.sqrt(jnp.diag(error_squared))
print(f"via autodiff:\nyield: {model.expected_data(result, include_auxdata=False)}\nerror: {error}")