Empirical Test Statistics

In this notebook we will compute test statistics empirically from pseudo-experiment and establish that they behave as assumed in the asymptotic approximation.

import numpy as np
import matplotlib.pyplot as plt
import pyhf

np.random.seed(0)
plt.rcParams.update({"font.size": 14})

First, we define the statistical model we will study.

• the signal expected event rate is 10 events

• the background expected event rate is 100 events

• a 10% uncertainty is assigned to the background

The expected event rates are chosen to lie comfortably in the asymptotic regime.

model = pyhf.simplemodels.uncorrelated_background(
    signal=[10.0], bkg=[100.0], bkg_uncertainty=[10.0]
)

The test statistics based on the profile likelihood described in arXiv:1007.1727 cover scanarios for both

• POI allowed to float to negative values (unbounded; \(\mu \in [-10, 10]\))

• POI constrained to non-negative values (bounded; \(\mu \in [0,10]\))

For consistency, test statistics (\(t_\mu, q_\mu\)) associated with bounded POIs are usually denoted with a tilde (\(\tilde{t}_\mu, \tilde{q}_\mu\)).

We set up the bounds for the fit as follows

unbounded_bounds = model.config.suggested_bounds()
unbounded_bounds[model.config.poi_index] = (-10, 10)

bounded_bounds = model.config.suggested_bounds()

Next we draw some synthetic datasets (also referrerd to as “toys” or pseudo-experiments). We will “throw” 300 toys:

true_poi = 1.0
n_toys = 300
toys = model.make_pdf(pyhf.tensorlib.astensor([true_poi, 1.0])).sample((n_toys,))

In the asymptotic treatment the test statistics are described as a function of the data’s best-fit POI value \(\hat\mu\).

So let’s run some fits so we can plots the empirical test statistics against \(\hat\mu\) to observed the emergence of the asymptotic behavior.

pars = np.asarray(
    [pyhf.infer.mle.fit(toy, model, par_bounds=unbounded_bounds) for toy in toys]
)
fixed_params = model.config.suggested_fixed()

We can now calculate all four test statistics described in arXiv:1007.1727

test_poi = 1.0
tmu = np.asarray(
    [
        pyhf.infer.test_statistics.tmu(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=unbounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

tmu_tilde = np.asarray(
    [
        pyhf.infer.test_statistics.tmu_tilde(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=bounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

qmu = np.asarray(
    [
        pyhf.infer.test_statistics.qmu(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=unbounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

qmu_tilde = np.asarray(
    [
        pyhf.infer.test_statistics.qmu_tilde(
            test_poi,
            toy,
            model,
            init_pars=model.config.suggested_init(),
            par_bounds=bounded_bounds,
            fixed_params=fixed_params,
        )
        for toy in toys
    ]
)

Let’s plot all the test statistics we have computed

muhat = pars[:, model.config.poi_index]
muhat_sigma = np.std(muhat)

We can check the asymptotic assumption that \(\hat{\mu}\) is distributed normally around it’s true value \(\mu' = 1\)

fig, ax = plt.subplots()
fig.set_size_inches(7, 5)

ax.set_xlabel(r"$\hat{\mu}$")
ax.set_ylabel("Density")
ax.set_ylim(top=0.5)

ax.hist(muhat, bins=np.linspace(-4, 5, 31), density=True)
ax.axvline(true_poi, label="true poi", color="black", linestyle="dashed")
ax.axvline(np.mean(muhat), label="empirical mean", color="red", linestyle="dashed")
ax.legend();

Here we define the asymptotic profile likelihood test statistics:

\[t_\mu = \frac{(\mu-\hat\mu)^2}{\sigma^2}\]

\[\begin{split}\tilde{t}_\mu = \begin{cases} t_\mu,\;\text{$\hat{\mu}>0$}\\ t_\mu - t_0,\; \text{else} \end{cases}\end{split}\]

\[\begin{split}q_\mu = \begin{cases} t_\mu,\;\text{$\hat{\mu}<\mu$}\\ 0,\; \text{else} \end{cases}\end{split}\]

\[\begin{split}\tilde{q}_\mu = \begin{cases} \tilde{t}_\mu,\;\text{$\hat{\mu}<\mu$}\\ 0,\; \text{else} \end{cases}\end{split}\]

def tmu_asymp(mutest, muhat, sigma):
    return (mutest - muhat) ** 2 / sigma**2


def tmu_tilde_asymp(mutest, muhat, sigma):
    a = tmu_asymp(mutest, muhat, sigma)
    b = tmu_asymp(mutest, muhat, sigma) - tmu_asymp(0.0, muhat, sigma)
    return np.where(muhat > 0, a, b)


def qmu_asymp(mutest, muhat, sigma):
    return np.where(
        muhat < mutest, tmu_asymp(mutest, muhat, sigma), np.zeros_like(muhat)
    )


def qmu_tilde_asymp(mutest, muhat, sigma):
    return np.where(
        muhat < mutest, tmu_tilde_asymp(mutest, muhat, sigma), np.zeros_like(muhat)
    )

And now we can compare them to the empirical values:

muhat_asymp = np.linspace(-3, 5)
fig, axarr = plt.subplots(2, 2)
fig.set_size_inches(14, 10)

labels = [r"$t_{\mu}$", "$\\tilde{t}_{\\mu}$", r"$q_{\mu}$", "$\\tilde{q}_{\\mu}$"]
data = [
    (tmu, tmu_asymp),
    (tmu_tilde, tmu_tilde_asymp),
    (qmu, qmu_asymp),
    (qmu_tilde, qmu_tilde_asymp),
]

for ax, (label, d) in zip(axarr.ravel(), zip(labels, data)):
    empirical, asymp_func = d
    ax.scatter(muhat, empirical, alpha=0.2, label=label)
    ax.plot(
        muhat_asymp,
        asymp_func(1.0, muhat_asymp, muhat_sigma),
        label=f"{label} asymptotic",
        c="r",
    )
    ax.set_xlabel(r"$\hat{\mu}$")
    ax.set_ylabel(f"{label}")
    ax.legend(loc="best")