Skip to content
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.

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

feat: Support configuring fixed values automatically from Model #1051

Merged
merged 24 commits into from
Sep 7, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
c31919b
fixed_vals add as argument to infer functions
cgottard Apr 26, 2020
80dafdb
fixed_vals add as argument to infer functions
cgottard Apr 26, 2020
32a7955
fixed_vals add as argument to infer functions
cgottard Apr 26, 2020
eb90ed5
tested
cgottard Apr 26, 2020
bf72f33
black formatting
cgottard Apr 26, 2020
5f8e16f
added fixed_vals positional argument to test backends
cgottard Apr 26, 2020
9437118
Adding NP fix and shift in test_backend_consistency
cgottard Apr 27, 2020
40b781b
correct histosys data structure in test_backend_consistency
cgottard Apr 27, 2020
bcd477d
fix test_backend_consistency
cgottard Apr 27, 2020
cc774e7
NP index fix in test_backend_consistency.py
cgottard May 5, 2020
d9e1972
Need 5 percent tolerane for 500 bin fit with TF
cgottard May 5, 2020
9d45958
pass in fixed_vals by converting to fixed_params
kratsg Sep 3, 2020
d3c3b47
make sure we propagate fixed_vals correctly everywhere
kratsg Sep 3, 2020
e0c82d1
fix up doc a bit
kratsg Sep 3, 2020
0c232a6
fix up notational difference
kratsg Sep 3, 2020
aada2e9
undo changes to tests for now
kratsg Sep 4, 2020
61b97ec
fixed it up
kratsg Sep 4, 2020
702b217
get tests passing
kratsg Sep 4, 2020
680502b
update a test and provide a fixed_params to deal with the pdf issue
kratsg Sep 4, 2020
d6831b6
scipy is able to minimize better so bump to 1 maxiter to force it to …
kratsg Sep 4, 2020
91847be
switch to fixed_params for public-facing calls in pyhf.infer
kratsg Sep 4, 2020
75fff60
update tests
kratsg Sep 4, 2020
a3cbe9f
fix up notebooks
kratsg Sep 4, 2020
3364a79
update fixed_poi_fit() to call fit() essentially
kratsg Sep 4, 2020
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
11 changes: 10 additions & 1 deletion docs/examples/notebooks/ImpactPlot.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,17 @@
" _, model, data = make_model([\"CRtt_meff\"])\n",
"\n",
" pyhf.set_backend(\"numpy\", pyhf.optimize.minuit_optimizer(verbose=True))\n",
" \n",
" constraints = constraints or []\n",
" init_pars = model.config.suggested_init()\n",
" fixed_params = model.config.suggested_fixed()\n",
" \n",
" for idx,fixed_val in constraints:\n",
" init_pars[idx] = fixed_val\n",
" fixed_params[idx] = True\n",
" \n",
" result = pyhf.infer.mle.fit(\n",
" data, model, fixed_vals=constraints, return_uncertainties=True\n",
" data, model, init_pars=init_pars, fixed_params=fixed_params, return_uncertainties=True\n",
" )\n",
" bestfit = result[:, 0]\n",
" errors = result[:, 1]\n",
Expand Down
11 changes: 10 additions & 1 deletion docs/examples/notebooks/pullplot.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -72,8 +72,17 @@
" _, model, data = make_model([\"CRtt_meff\"])\n",
"\n",
" pyhf.set_backend(\"numpy\", pyhf.optimize.minuit_optimizer(verbose=True))\n",
" \n",
" constraints = constraints or []\n",
" init_pars = model.config.suggested_init()\n",
" fixed_params = model.config.suggested_fixed()\n",
" \n",
" for idx,fixed_val in constraints:\n",
" init_pars[idx] = fixed_val\n",
" fixed_params[idx] = True\n",
" \n",
" result = pyhf.infer.mle.fit(\n",
" data, model, fixed_vals=constraints, return_uncertainties=True\n",
" data, model, init_pars=init_pars, fixed_params=fixed_params, return_uncertainties=True\n",
" )\n",
" bestfit = result[:, 0]\n",
" errors = result[:, 1]\n",
Expand Down
22 changes: 17 additions & 5 deletions src/pyhf/infer/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,14 @@


def hypotest(
poi_test, data, pdf, init_pars=None, par_bounds=None, qtilde=False, **kwargs
poi_test,
data,
pdf,
init_pars=None,
par_bounds=None,
fixed_params=None,
qtilde=False,
**kwargs,
):
r"""
Compute :math:`p`-values and test statistics for a single value of the parameter of interest.
Expand All @@ -32,8 +39,9 @@ def hypotest(
poi_test (Number or Tensor): The value of the parameter of interest (POI)
data (Number or Tensor): The data considered
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``
init_pars (Array or Tensor): The initial parameter values to be used for minimization
par_bounds (Array or Tensor): The parameter value bounds to be used for minimization
init_pars (`tensor`): The initial parameter values to be used for minimization
par_bounds (`tensor`): The parameter value bounds to be used for minimization
fixed_params (`tensor`): Whether to fix the parameter to the init_pars value during minimization
qtilde (Bool): When ``True`` perform the calculation using the alternative
test statistic, :math:`\tilde{q}_{\mu}`, as defined under the Wald
approximation in Equation (62) of :xref:`arXiv:1007.1727`.
Expand Down Expand Up @@ -118,15 +126,19 @@ def hypotest(
"""
init_pars = init_pars or pdf.config.suggested_init()
par_bounds = par_bounds or pdf.config.suggested_bounds()
tensorlib, _ = get_backend()
fixed_params = fixed_params or pdf.config.suggested_fixed()

calc = AsymptoticCalculator(data, pdf, init_pars, par_bounds, qtilde=qtilde)
calc = AsymptoticCalculator(
data, pdf, init_pars, par_bounds, fixed_params, qtilde=qtilde
)
teststat = calc.teststatistic(poi_test)
sig_plus_bkg_distribution, b_only_distribution = calc.distributions(poi_test)

CLsb = sig_plus_bkg_distribution.pvalue(teststat)
CLb = b_only_distribution.pvalue(teststat)
CLs = CLsb / CLb

tensorlib, _ = get_backend()
# Ensure that all CL values are 0-d tensors
CLsb, CLb, CLs = (
tensorlib.astensor(CLsb),
Expand Down
41 changes: 35 additions & 6 deletions src/pyhf/infer/calculators.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from .test_statistics import qmu, qmu_tilde


def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_params):
"""
Compute Asimov Dataset (expected yields at best-fit values) for a given POI value.

Expand All @@ -22,12 +22,15 @@ def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
init_pars (`tensor`): The initial parameter values to be used for fitting.
par_bounds (`tensor`): The parameter value bounds to be used for fitting.
fixed_params (`tensor`): Parameters to be held constant in the fit.

Returns:
Tensor: The Asimov dataset.

"""
bestfit_nuisance_asimov = fixed_poi_fit(asimov_mu, data, pdf, init_pars, par_bounds)
bestfit_nuisance_asimov = fixed_poi_fit(
asimov_mu, data, pdf, init_pars, par_bounds, fixed_params
)
return pdf.expected_data(bestfit_nuisance_asimov)


Expand Down Expand Up @@ -118,7 +121,15 @@ def expected_value(self, nsigma):
class AsymptoticCalculator(object):
"""The Asymptotic Calculator."""

def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
def __init__(
self,
data,
pdf,
init_pars=None,
par_bounds=None,
fixed_params=None,
qtilde=False,
):
"""
Asymptotic Calculator.

Expand All @@ -127,6 +138,7 @@ def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
init_pars (`tensor`): The initial parameter values to be used for fitting.
par_bounds (`tensor`): The parameter value bounds to be used for fitting.
fixed_params (`tensor`): Whether to fix the parameter to the init_pars value during minimization
qtilde (`bool`): Whether to use qtilde as the test statistic.

Returns:
Expand All @@ -137,6 +149,8 @@ def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
self.pdf = pdf
self.init_pars = init_pars or pdf.config.suggested_init()
self.par_bounds = par_bounds or pdf.config.suggested_bounds()
self.fixed_params = fixed_params or pdf.config.suggested_fixed()

self.qtilde = qtilde
self.sqrtqmuA_v = None

Expand Down Expand Up @@ -173,16 +187,31 @@ def teststatistic(self, poi_test):
teststat_func = qmu_tilde if self.qtilde else qmu

qmu_v = teststat_func(
poi_test, self.data, self.pdf, self.init_pars, self.par_bounds
poi_test,
self.data,
self.pdf,
self.init_pars,
self.par_bounds,
self.fixed_params,
)
sqrtqmu_v = tensorlib.sqrt(qmu_v)

asimov_mu = 0.0
asimov_data = generate_asimov_data(
asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds
asimov_mu,
self.data,
self.pdf,
self.init_pars,
self.par_bounds,
self.fixed_params,
)
qmuA_v = teststat_func(
poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds
poi_test,
asimov_data,
self.pdf,
self.init_pars,
self.par_bounds,
self.fixed_params,
)
self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v)

Expand Down
43 changes: 27 additions & 16 deletions src/pyhf/infer/mle.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,9 @@ def twice_nll(pars, data, pdf):
return -2 * pdf.logpdf(pars, data)


def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
def fit(data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs):
r"""
Run a unconstrained maximum likelihood fit.
Run a maximum likelihood fit.
This is done by minimizing the objective function :func:`~pyhf.infer.mle.twice_nll`
of the model parameters given the observed data.
This is used to produce the maximal likelihood :math:`L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)`
Expand Down Expand Up @@ -87,6 +87,7 @@ def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
init_pars (`list`): Values to initialize the model parameters at for the fit
par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
fixed_params (`list`): Parameters to be held constant in the fit.
kwargs: Keyword arguments passed through to the optimizer API

Returns:
Expand All @@ -96,10 +97,23 @@ def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
_, opt = get_backend()
init_pars = init_pars or pdf.config.suggested_init()
par_bounds = par_bounds or pdf.config.suggested_bounds()
return opt.minimize(twice_nll, data, pdf, init_pars, par_bounds, **kwargs)
fixed_params = fixed_params or pdf.config.suggested_fixed()

# get fixed vals from the model
fixed_vals = [
(index, init)
for index, (init, is_fixed) in enumerate(zip(init_pars, fixed_params))
if is_fixed
]

def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs):
return opt.minimize(
twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs
)


def fixed_poi_fit(
poi_val, data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs
):
r"""
Run a maximum likelihood fit with the POI value fixed.
This is done by minimizing the objective function of :func:`~pyhf.infer.mle.twice_nll`
Expand Down Expand Up @@ -142,6 +156,7 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs)
pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
init_pars (`list`): Values to initialize the model parameters at for the fit
par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
fixed_params (`list`): Parameters to be held constant in the fit.
kwargs: Keyword arguments passed through to the optimizer API

Returns:
Expand All @@ -152,15 +167,11 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs)
raise UnspecifiedPOI(
'No POI is defined. A POI is required to fit with a fixed POI.'
)
_, opt = get_backend()
init_pars = init_pars or pdf.config.suggested_init()
par_bounds = par_bounds or pdf.config.suggested_bounds()
return opt.minimize(
twice_nll,
data,
pdf,
init_pars,
par_bounds,
[(pdf.config.poi_index, poi_val)],
**kwargs,
)

init_pars = [*(init_pars or pdf.config.suggested_init())]
fixed_params = [*(fixed_params or pdf.config.suggested_fixed())]

init_pars[pdf.config.poi_index] = poi_val
fixed_params[pdf.config.poi_index] = True

return fit(data, pdf, init_pars, par_bounds, fixed_params, **kwargs)
Loading