refactor: Use MLE API (#690)

lukasheinrich · matthewfeickert · commit 51d8a39086ff · 2019-12-24T14:26:32.000-06:00
* Simplify optimizer API to only have 'minimize()'
* Add option to return minimized function in 'minimize()'
* Add option to return uncertainties in the parameters for 'minimize()'
* Move concrete fits that are needed for maximum likelihood estimation / inference to pyhf.infer.mle
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -99,7 +99,7 @@ jobs:
     - name: Check docstrings
       run: |
         # Group 1 is related to docstrings
-        pydocstyle --select D1 src/pyhf/pdf.py src/pyhf/probability.py src/pyhf/interpolators
+        pydocstyle --select D1 src/pyhf/pdf.py src/pyhf/probability.py src/pyhf/interpolators src/pyhf/infer src/pyhf/optimize
     - name: Test and build docs
       run: |
         python -m doctest README.md
diff --git a/docs/api.rst b/docs/api.rst
@@ -119,7 +119,9 @@ Inference
 
    hypotest
    test_statistics.qmu
-   utils.loglambdav
+   mle.twice_nll
+   mle.fit
+   mle.fixed_poi_fit
    utils.generate_asimov_data
    utils.pvals_from_teststat
    utils.pvals_from_teststat_expected
diff --git a/docs/examples/notebooks/ImpactPlot.ipynb b/docs/examples/notebooks/ImpactPlot.ipynb
diff --git a/docs/examples/notebooks/ShapeFactor.ipynb b/docs/examples/notebooks/ShapeFactor.ipynb
@@ -181,8 +181,7 @@
    "source": [
     "print('initialization parameters: {}'.format(pdf.config.suggested_init()))\n",
     "\n",
-    "unconpars = pyhf.optimizer.unconstrained_bestfit(pyhf.infer.utils.loglambdav, data, pdf,\n",
-    "                                                 pdf.config.suggested_init(), pdf.config.suggested_bounds())\n",
+    "unconpars = pyhf.infer.mle.fit(data, pdf)\n",
     "print('parameters post unconstrained fit: {}'.format(unconpars))"
    ]
   },
@@ -284,4 +283,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
diff --git a/docs/examples/notebooks/binderexample/StatisticalAnalysis.ipynb b/docs/examples/notebooks/binderexample/StatisticalAnalysis.ipynb
@@ -1072,8 +1072,7 @@
     "nominal = pdf.config.suggested_init()\n",
     "background_only = pdf.config.suggested_init()\n",
     "background_only[pdf.config.poi_index] = 0.0\n",
-    "best_fit = pyhf.optimizer.unconstrained_bestfit(\n",
-    "    pyhf.infer.utils.loglambdav, data, pdf, pdf.config.suggested_init(), pdf.config.suggested_bounds())"
+    "best_fit = pyhf.infer.mle.fit(data, pdf)"
    ]
   },
   {
@@ -9051,4 +9050,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}
diff --git a/docs/examples/notebooks/multiBinPois.ipynb b/docs/examples/notebooks/multiBinPois.ipynb
@@ -200,7 +200,7 @@
     "\n",
     "print(init_pars)\n",
     "\n",
-    "bestfit_pars = optimizer.unconstrained_bestfit(pyhf.infer.utils.loglambdav, data, pdf, init_pars, par_bounds)\n",
+    "bestfit_pars = pyhf.infer.mle.fit(data, pdf, init_pars, par_bounds)\n",
     "bestfit_cts  = pdf.expected_data(bestfit_pars, include_auxdata = False)"
    ]
   },
@@ -371,4 +371,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 1
-}
+}
diff --git a/docs/examples/notebooks/multichannel-coupled-histo.ipynb b/docs/examples/notebooks/multichannel-coupled-histo.ipynb
diff --git a/docs/examples/notebooks/pullplot.ipynb b/docs/examples/notebooks/pullplot.ipynb
diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py
@@ -0,0 +1,63 @@
+"""Module for Maximum Likelihood Estimation."""
+from .. import get_backend
+
+
+def twice_nll(pars, data, pdf):
+    """
+    Twice the negative Log-Likelihood.
+
+    Args:
+        data (`tensor`): the data
+        pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
+
+    Returns:
+        Twice the negative log likelihood.
+
+    """
+    return -2 * pdf.logpdf(pars, data)
+
+
+def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
+    """
+    Run a unconstrained maximum likelihood fit.
+
+    Args:
+        data (`tensor`): the data
+        pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
+        kwargs: keyword arguments passed through to the optimizer API
+
+    Returns:
+        see optimizer API
+
+    """
+    _, 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)
+
+
+def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs):
+    """
+    Run a maximum likelihood fit with the POI value fixzed.
+
+    Args:
+        data: the data
+        pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
+        kwargs: keyword arguments passed through to the optimizer API
+
+    Returns:
+        see optimizer API
+
+    """
+    _, 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
+    )
diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py
@@ -1,5 +1,5 @@
 from .. import get_backend
-from .utils import loglambdav
+from .mle import fixed_poi_fit, fit
 
 
 def qmu(mu, data, pdf, init_pars, par_bounds):
@@ -30,13 +30,13 @@ def qmu(mu, data, pdf, init_pars, par_bounds):
         Float: The calculated test statistic, :math:`q_{\mu}`
     """
     tensorlib, optimizer = get_backend()
-    mubhathat = optimizer.constrained_bestfit(
-        loglambdav, mu, data, pdf, init_pars, par_bounds
+    mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit(
+        mu, data, pdf, init_pars, par_bounds, return_fitted_val=True
     )
-    muhatbhat = optimizer.unconstrained_bestfit(
-        loglambdav, data, pdf, init_pars, par_bounds
+    muhatbhat, unconstrained_fit_lhood_val = fit(
+        data, pdf, init_pars, par_bounds, return_fitted_val=True
     )
-    qmu = loglambdav(mubhathat, data, pdf) - loglambdav(muhatbhat, data, pdf)
+    qmu = fixed_poi_fit_lhood_val - unconstrained_fit_lhood_val
     qmu = tensorlib.where(
         muhatbhat[pdf.config.poi_index] > mu, tensorlib.astensor([0]), qmu
     )
diff --git a/src/pyhf/infer/utils.py b/src/pyhf/infer/utils.py
@@ -1,16 +1,11 @@
 """Utility Functions for model inference."""
 from .. import get_backend
-
-
-def loglambdav(pars, data, pdf):
-    return -2 * pdf.logpdf(pars, data)
+from .mle import fixed_poi_fit
 
 
 def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
-    _, optimizer = get_backend()
-    bestfit_nuisance_asimov = optimizer.constrained_bestfit(
-        loglambdav, asimov_mu, data, pdf, init_pars, par_bounds
-    )
+    """Compute Asimov Dataset (expected yields at best-fit values) for a given POI value."""
+    bestfit_nuisance_asimov = fixed_poi_fit(asimov_mu, data, pdf, init_pars, par_bounds)
     return pdf.expected_data(bestfit_nuisance_asimov)
 
 
diff --git a/src/pyhf/optimize/__init__.py b/src/pyhf/optimize/__init__.py
@@ -1,3 +1,5 @@
+"""Optimizers for Tensor Functions."""
+
 from .. import exceptions
 
 
diff --git a/src/pyhf/optimize/autodiff.py b/src/pyhf/optimize/autodiff.py
@@ -1,9 +1,31 @@
+"""Helper Classes for use of automatic differentiation."""
 import scipy
 from .. import get_backend
 
 
 class AutoDiffOptimizerMixin(object):
-    def minimize(self, objective, data, pdf, init_pars, par_bounds, fixed_vals=None):
+    """Mixin Class to build optimizers that use automatic differentiation."""
+
+    def __init__(*args, **kwargs):
+        """Create Mixin for autodiff-based optimizers."""
+
+    def minimize(
+        self,
+        objective,
+        data,
+        pdf,
+        init_pars,
+        par_bounds,
+        fixed_vals=None,
+        return_fitted_val=False,
+    ):
+        """
+        Find Function Parameters that minimize the Objective.
+
+        Returns:
+            bestfit parameters
+        
+        """
         tensorlib, _ = get_backend()
         tv, fixed_values_tensor, func, init, bounds = self.setup_minimize(
             objective, data, pdf, init_pars, par_bounds, fixed_vals
@@ -12,19 +34,10 @@ def minimize(self, objective, data, pdf, init_pars, par_bounds, fixed_vals=None)
             func, init, method='SLSQP', jac=True, bounds=bounds
         )
         nonfixed_vals = fitresult.x
-        return tv.stitch([fixed_values_tensor, tensorlib.astensor(nonfixed_vals)])
-
-    def unconstrained_bestfit(self, objective, data, pdf, init_pars, par_bounds):
-        return self.minimize(objective, data, pdf, init_pars, par_bounds)
-
-    def constrained_bestfit(
-        self, objective, constrained_mu, data, pdf, init_pars, par_bounds
-    ):
-        return self.minimize(
-            objective,
-            data,
-            pdf,
-            init_pars,
-            par_bounds,
-            [(pdf.config.poi_index, constrained_mu)],
+        fitted_val = fitresult.fun
+        fitted_pars = tv.stitch(
+            [fixed_values_tensor, tensorlib.astensor(nonfixed_vals)]
         )
+        if return_fitted_val:
+            return fitted_pars, tensorlib.astensor(fitted_val)
+        return fitted_pars
diff --git a/src/pyhf/optimize/opt_minuit.py b/src/pyhf/optimize/opt_minuit.py
@@ -1,3 +1,5 @@
+"""MINUIT Optimizer Backend."""
+
 import iminuit
 import logging
 import numpy as np
@@ -6,7 +8,16 @@
 
 
 class minuit_optimizer(object):
+    """MINUIT Optimizer Backend."""
+
     def __init__(self, verbose=False, ncall=10000, errordef=1, steps=1000):
+        """
+        Create MINUIT Optimizer.
+
+        Args:
+            verbose (`bool`): print verbose output during minimization
+        
+        """
         self.verbose = 0
         self.ncall = ncall
         self.errordef = errordef
@@ -45,23 +56,32 @@ def f(pars):
         )
         return mm
 
-    def minimize(self, objective, data, pdf, init_pars, par_bounds, fixed_vals=None):
+    def minimize(
+        self,
+        objective,
+        data,
+        pdf,
+        init_pars,
+        par_bounds,
+        fixed_vals=None,
+        return_fitted_val=False,
+        return_uncertainties=False,
+    ):
+        """
+        Find Function Parameters that minimize the Objective.
+
+        Returns:
+            bestfit parameters
+        
+        """
         mm = self._make_minuit(objective, data, pdf, init_pars, par_bounds, fixed_vals)
         result = mm.migrad(ncall=self.ncall)
         assert result
-        return np.asarray([x[1] for x in mm.values.items()])
-
-    def unconstrained_bestfit(self, objective, data, pdf, init_pars, par_bounds):
-        return self.minimize(objective, data, pdf, init_pars, par_bounds)
-
-    def constrained_bestfit(
-        self, objective, constrained_mu, data, pdf, init_pars, par_bounds
-    ):
-        return self.minimize(
-            objective,
-            data,
-            pdf,
-            init_pars,
-            par_bounds,
-            [(pdf.config.poi_index, constrained_mu)],
-        )
+        if return_uncertainties:
+            bestfit_pars = np.asarray([(v, mm.errors[k]) for k, v in mm.values.items()])
+        else:
+            bestfit_pars = np.asarray([v for k, v in mm.values.items()])
+        bestfit_value = mm.fval
+        if return_fitted_val:
+            return bestfit_pars, bestfit_value
+        return bestfit_pars
diff --git a/src/pyhf/optimize/opt_pytorch.py b/src/pyhf/optimize/opt_pytorch.py
@@ -1,16 +1,29 @@
+"""PyTorch Optimizer Backend."""
+
 from .. import get_backend, default_backend
 from ..tensor.common import _TensorViewer
 from .autodiff import AutoDiffOptimizerMixin
 import torch
 
 
 class pytorch_optimizer(AutoDiffOptimizerMixin):
-    def __init__(self, *args, **kargs):
-        pass
+    """PyTorch Optimizer Backend."""
 
     def setup_minimize(
         self, objective, data, pdf, init_pars, par_bounds, fixed_vals=None
     ):
+        """
+        Prepare Minimization for AutoDiff-Optimizer.
+
+        Args:
+            objective: objective function
+            data: observed data
+            pdf: model
+            init_pars: initial parameters
+            par_bounds: parameter boundaries
+            fixed_vals: fixed parameter values
+
+        """
         tensorlib, _ = get_backend()
         all_idx = default_backend.astensor(range(pdf.config.npars), dtype='int')
         all_init = default_backend.astensor(init_pars)
diff --git a/src/pyhf/optimize/opt_scipy.py b/src/pyhf/optimize/opt_scipy.py
diff --git a/src/pyhf/optimize/opt_tflow.py b/src/pyhf/optimize/opt_tflow.py
diff --git a/tests/test_optim.py b/tests/test_optim.py
