Upper limit calculation

[gollumben]

Dear pyhf team,

thanks again for providing this wonderful software suite with the marvelous documentation.

In the standard pyhf implementation, the upper limit is found by scanning over a range of mu_sig and interpolating between the tested points. Since it takes ~3 minutes to calculate one CLs value for the workspace under consideration, this approach is time consuming, if running with high granularity. If running with low granularity, the result is rather different compared to the limits obtained with this asymptotic calculator tool (please let me know if you do not have access) (11 points between mu=0 and mu=5 scanned):

• Asymptotic tool: mu_lim=0.40
• pyhf: mu_lim=0.49

The time consumption of the asymptotic calculator and pyhf is comparable (~30 minutes) on my machine for the attached workspace. These values might seem reasonably close, but for other points the observed difference is larger (e.g. 0.2 VS 0.4; I suppose it is due to the coarse segmentation). Furthermore, I get these warnings, which I suppose come from NaNs.

RuntimeWarning: invalid value encountered in greater_equal

I also tested using the autoscan method implemented in #1274. Unfortunately, I get a NaN error for the workspace attached.

Do you have an idea what I am doing wrong or suggestions to speed up the calculations? In the end I will want to run on multiple points with mu_lim=0.1-200, so a wide range would be useful. Could it be beneficial to port the asymptotic tool to pyhf (I think a lot of things in the script are already in pyhf or not needed)?

I am attaching the code used to calculate the limits below. This discussion could be related to #1695.

Looking forward to you feedback and thanks in advance,
Ben

---

Attachments:

0L1LB_a100_DM10_H1000_tb1_st0p7.json.txt

import json, pyhf
import numpy as np

def runExclFit(workspaceJsonFileName, mode):
	print("Running excl. fit on:",workspaceJsonFileName,"with mode:",mode)
	pyhf.set_backend("numpy")

	# open workspace
	with open(workspaceJsonFileName) as iF:
		jsonObj = json.load(iF)

	workspace = pyhf.Workspace(jsonObj)

	model = workspace.model(measurement_name="NormalMeasurement")
	data = workspace.data(model)

	# calculate CLs
	if mode == "CLs":
		mu_test = 1.0 ## threshold for CLs
		CLs_obs, CLs_tails, CLs_exp_band = pyhf.infer.hypotest(
			mu_test, data, model, return_tail_probs = True, return_expected_set=True, test_stat="qtilde"
		)
		print("CLs_obs, CLs_tails, CLs_exp_band",CLs_obs,CLs_tails,CLs_exp_band)

	# calculate upper limit
	elif "Upperlimit" in mode:
		if "lin" in mode:
			scan = np.linspace(0, 5, 11) ## to be used with fixed-range mu scan range
		else:
			scan = "auto" ## automatic range determination
		obs_limit, exp_limits = pyhf.infer.intervals.upperlimit(
			data, model, scan
		)
		print("obs_limit, exp_limits",obs_limit,exp_limits)
[0L1LB_a100_DM10_H1000_tb1_st0p7.json.txt](https://github.com/scikit-hep/pyhf/files/7649001/0L1LB_a100_DM10_H1000_tb1_st0p7.json.txt)


if __name__ == "__main__":
	filename = "/lustre/fs22/group/atlas/brueers/tW0LBoosted/combination2L/noCorrNP/0L1LB_a100_DM10_H1000_tb1_st0p7.json"
	
	## value from HF: 0.00
	runExclFit(filename, "CLs") # gives 1.95e-5 --> compatible

	## value from "Asymptotic Calculator" script: 0.40
	runExclFit(filename, "Upperlimit_lin") # gives 0.47
	runExclFit(filename, "Upperlimit") # gives nan and fails 

[alexander-held]

The script you linked uses some clever tricks to find the next point and converge very efficiently. If speed is crucial, those probably would also need to be included. I think the currently available grid scan was never meant to be efficient and cannot compete with highly optimized implementations. Even the autoscan version may need some further tuning to converge in as few steps as that script.

Regarding the NaNs, I think there can be (at least) two scenarios where they appear:

• workspace predicts negative yields for some parameter settings, possibly causing fits to fail (see Fit succeeds with scipy, but fails with minuit #1695 (comment), I ran a simple MLE fit on the workspace provided here and again see NaN likelihoods for some regions of parameter space),
• the initial POI value tested is very small, resulting in CLs = NaN (that could be an issue if you expect a wide range of possible limits).