Skip to content

Implement observed asymptotic CLs likelihood limit #136

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
davidwalter2 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Implement observed asymptotic CLs likelihood limit #136

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 58 additions & 1 deletion .github/workflows/main.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -429,11 +429,68 @@ jobs:
rabbit_limit.py $RABBIT_OUTDIR/test_tensor_bsm.hdf5 -o $RABBIT_OUTDIR/ --postfix bsm
-t 0 --unblind --expectSignal bsm1 0 --expectSignal bsm2 1 --allowNegativeParam --asymptoticLimits bsm1 --levels 0.05

- name: bsm limits on channel
- name: check bsm observed likelihood limit
run: |
python - <<'PY'
import os
import numpy as np
from rabbit import io_tools

fr = io_tools.get_fitresult(
os.path.join(os.environ["RABBIT_OUTDIR"], "limits_bsm.hdf5")
)
keys = list(fr.keys())
assert "likelihood_asymptoticLimits_observed" in keys, keys
obs = np.array(
fr["likelihood_asymptoticLimits_observed"].get().values()
).flatten()
assert np.all(np.isfinite(obs)) and np.all(obs > 0), obs
# observed must sit inside the expected +/-2 sigma band (clb axis
# ordered 0.025 .. 0.975), with padding for the gaussian/likelihood
# method difference; this catches NaN / sign / scale regressions
exp = np.array(fr["likelihood_asymptoticLimits_expected"].get().values())
lo = exp[..., 0].flatten()
hi = exp[..., -1].flatten()
assert np.all(obs >= 0.5 * lo) and np.all(obs <= 2.0 * hi), (obs, lo, hi)
gobs = np.array(
fr["gaussian_asymptoticLimits_observed"].get().values()
).flatten()
assert np.all(obs > 0.3 * gobs) and np.all(obs < 3.0 * gobs), (obs, gobs)
print("likelihood observed limit OK:", obs)
PY

- name: bsm limits on channel
run: >-
rabbit_limit.py $RABBIT_OUTDIR/test_tensor_bsm.hdf5 -o $RABBIT_OUTDIR/ --postfix bsm
-t 0 --unblind --paramModel Mixture bsm1 sig --expectSignal bsm1_sig_mixing 0 --asymptoticLimits bsm1_sig_mixing --levels 0.05 --allowNegativeParam

- name: check bsm observed likelihood limit on channel
run: |
python - <<'PY'
import os
import numpy as np
from rabbit import io_tools

fr = io_tools.get_fitresult(
os.path.join(os.environ["RABBIT_OUTDIR"], "limits_bsm.hdf5")
)
keys = list(fr.keys())
assert "likelihood_asymptoticLimits_observed" in keys, keys
obs = np.array(
fr["likelihood_asymptoticLimits_observed"].get().values()
).flatten()
assert np.all(np.isfinite(obs)) and np.all(obs > 0), obs
exp = np.array(fr["likelihood_asymptoticLimits_expected"].get().values())
lo = exp[..., 0].flatten()
hi = exp[..., -1].flatten()
assert np.all(obs >= 0.5 * lo) and np.all(obs <= 2.0 * hi), (obs, lo, hi)
gobs = np.array(
fr["gaussian_asymptoticLimits_observed"].get().values()
).flatten()
assert np.all(obs > 0.3 * gobs) and np.all(obs < 3.0 * gobs), (obs, gobs)
print("likelihood observed limit on channel OK:", obs)
PY

plotting:
runs-on: [self-hosted, linux, x64]
needs: [setenv, fitting]
Expand Down
45 changes: 27 additions & 18 deletions bin/rabbit_limit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,6 +166,9 @@ def do_asymptotic_limits(args, fitter, ws, data_values, mapping=None, fit_data=F

# this is the denominator of q for likelihood based limits
nllvalreduced_asimov = fitter_asimov.reduced_nll().numpy()
if fit_data:
# global minimum on data (unconditional fit, done above)
nllvalreduced = fitter.reduced_nll().numpy()

for j, cl_s in enumerate(args.levels):
logger.info(f" -- AsymptoticLimits ( CLs={round(cl_s*100,1):4.1f}% ) -- ")
Expand Down Expand Up @@ -200,16 +203,24 @@ def do_asymptotic_limits(args, fitter, ws, data_values, mapping=None, fit_data=F
key, xobs, xobs_err, xerr, cl_s
)

if "likelihood" in args.modes:
# Likelihood approximation
if mapping is not None:
# TODO: implement for channels
pass
else:
# TODO: make it work
# nllvalreduced = fitter.reduced_nll().numpy()
# limits_nll_obs[i, j] = asymptotic_limits.compute_likelihood_limit(fitter, fitter_asimov, nllvalreduced, nllvalreduced_asimov, key, cl_s)
pass
if "likelihood" in args.modes and fit_data:
# Likelihood approximation (covers both the parameter case
# fun=None and the channel/mapping case fun=mapping.compute_flat)
limits_nll_obs[i, j] = asymptotic_limits.compute_likelihood_limit(
fitter,
fitter_asimov,
nllvalreduced,
nllvalreduced_asimov,
key,
cl_s,
fun=fun,
r_init=(
limits_obs[i, j]
if (mapping is None and "gaussian" in args.modes)
else None
),
allow_negative=args.allowNegativeParam,
)

if "gaussian" in args.modes:
ws.add_limits_hist(
Expand All @@ -228,7 +239,6 @@ def do_asymptotic_limits(args, fitter, ws, data_values, mapping=None, fit_data=F
)

if "likelihood" in args.modes:
# TODO: implement for channels
ws.add_limits_hist(
limits_nll,
args.asymptoticLimits,
Expand All @@ -237,13 +247,12 @@ def do_asymptotic_limits(args, fitter, ws, data_values, mapping=None, fit_data=F
base_name="likelihood_asymptoticLimits_expected",
)

# TODO: make it work
# ws.add_limits_hist(
# limits_nll_obs,
# args.asymptoticLimits,
# args.levels,
# base_name="likelihood_asymptoticLimits_observed",
# )
ws.add_limits_hist(
limits_nll_obs,
args.asymptoticLimits,
args.levels,
base_name="likelihood_asymptoticLimits_observed",
)


def main():
Expand Down
Loading