Can't find the likelihood maximum

[APMDSLHC]

System Settings

Fedora Linux 35 Python 3.9.12

Describe the bug

Optimizer in combiner.maximize_likelihood() can't reach the required precision.

Sometimes it crashes, sometimes it doesn't.

fun: nan
     jac: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan])
 message: 'Inequality constraints incompatible'
    nfev: 61
     nit: 1
    njev: 1
  status: 4
 success: False
       x: array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  1.        , -8.87822344,  1.        ,
        0.        ,  0.        ,  1.        ,  1.        ,  1.        ])
Traceback (most recent call last):
  File "/home/pascal/anaconda3/lib/python3.9/site-packages/pyhf/optimize/mixins.py", line 64, in _internal_minimize
    assert result.success
AssertionError
     fun: nan
     jac: array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
       nan, nan, nan, nan, nan, nan, nan, nan])
 message: 'Inequality constraints incompatible'
    nfev: 61
     nit: 1
    njev: 1
  status: 4
 success: False
       x: array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  1.        , -7.87822344,  1.        ,
        0.        ,  0.        ,  1.        ,  1.        ,  1.        ])
Traceback (most recent call last):
  File "/home/pascal/anaconda3/lib/python3.9/site-packages/pyhf/optimize/mixins.py", line 64, in _internal_minimize
    assert result.success
AssertionError

---------------------------------------------------------------------------
RuntimeWarning                            Traceback (most recent call last)
Input In [8], in <cell line: 7>()
      4 nll_observed = np.array([combiner.likelihood(mu=m, return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.observed) for m in mu])
      6 muhat_obs, nll_min_obs = combiner.maximize_likelihood(return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.observed)
----> 7 muhat_apri, nll_min_apri = combiner.maximize_likelihood(return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.apriori)
      9 fig = plt.figure(figsize=(8,6))
     11 plt.plot(mu, nll_apriori, label=str(madstats.ExpectationType.apriori), color="tab:red")

File ~/SModelS/madstats/src/madstats/combiner/prediction_combiner.py:156, in PredictionCombiner.maximize_likelihood(self, return_nll, expected, allow_negative_signal, iteration_threshold, **kwargs)
    147 opt = scipy.optimize.minimize(
    148     negloglikelihood,
    149     muhat_init,
   (...)
    152     options={"maxiter": iteration_threshold},
    153 )
    155 if not opt.success:
--> 156     raise RuntimeWarning("Optimiser was not able to reach required precision.")
    158 nll, muhat = opt.fun, opt.x[0]
    159 if not allow_negative_signal and muhat < 0.0:

RuntimeWarning: Optimiser was not able to reach required precision.

To Reproduce

import madstats, json
import numpy as np

with open("./madstats/test/signal_test.json", "r") as f:
    signal = json.load(f)

with open("./madstats/test/background_test.json", "r") as f:
    background = json.load(f)

stat_model = madstats.get_multi_region_statistical_model(
    analysis="atlas_susy_2018_31",
    signal=signal,
    background=background,
    xsection=0.000207244,
)

dat = np.load("./madstats/test/simplified_likelihood_data.npz")
sl_model = madstats.get_multi_region_statistical_model(
    analysis="cms_sus_19_006",
    signal=dat["nsignal"],
    background=dat["observed"],
    covariance=dat["covariance"],
    nb=dat["backgrounds"],
    xsection=0.000207244,
    delta_sys=0.
)

combiner = madstats.PredictionCombiner(sl_model, stat_model)

muhat_apri, nll_min_apri = combiner.maximize_likelihood(return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.apriori)

Expected behaviour

Finding the parameters that maximise the likelihood.

Additional information

numpy version 1.21.6 scipy version 1.8.0

[jackaraz]

Hi @APMDSLHC, I modified the code but left the warnings on for now. So it should warn you still, but you should still get a result. If you like, you can change the default backend-specific settings in the combiner e.g.

import madstats
kwargs = {str(madstats.AvailableBackends.pyhf): {"iteration_threshold": 20}}
muhat_apri, nll_min_apri = combiner.maximize_likelihood(
    return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.apriori, **kwargs
)

You can check the appropriate function's keywords to see what you can modify via kwargs e.g.

madstats.backends.pyhf_backend.interface.PyhfInterface.likelihood?

# Out:
Signature:
madstats.backends.pyhf_backend.interface.PyhfInterface.likelihood(
    self,
    mu: Optional[float] = 1.0,
    expected: Optional[madstats.utils.ExpectationType] = observed,
    allow_negative_signal: bool = True,
    return_nll: Optional[bool] = False,
    isAsimov: Optional[bool] = False,
    iteration_threshold: Optional[int] = 10,
) -> float
Docstring:
Compute the likelihood of the given statistical model

:param mu: POI (signal strength)
:param expected: observed, expected (true, apriori) or aposteriori
:param allow_negative_signal: if true, POI can get negative values
:param return_nll: if true returns negative log-likelihood value
:param isAsimov: if true, computes likelihood for Asimov data
:param mu_lim: boundaries for mu
:param iteration_threshold: number of iterations to be held for convergence of the fit.
:return: (float) likelihood
Type:      function

[APMDSLHC]

Thanks, but I now have infinities for some values of mu.

When I do:

with open("./madstats/test/signal_test.json", "r") as f:
    signal = json.load(f)

with open("./madstats/test/background_test.json", "r") as f:
    background = json.load(f)

stat_model = madstats.get_multi_region_statistical_model(
    analysis="atlas_susy_2018_31",
    signal=signal,
    background=background,
    xsection=0.000207244,
)

mu = np.linspace(-3,2,25)

nll_apriori = np.array([stat_model.likelihood(mu=m, return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.apriori) for m in mu])
nll_observed = np.array([stat_model.likelihood(mu=m, return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.observed) for m in mu])

muhat_obs, nll_min_obs = stat_model.maximize_likelihood(return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.observed)
muhat_apri, nll_min_apri = stat_model.maximize_likelihood(return_nll=True, allow_negative_signal=True, expected=madstats.ExpectationType.apriori)

print(nll_apriori)
print(nll_observed)

I get:

[           inf   358.90796301            inf  2007.40776346
            inf  8467.43623774            inf 22949.41643341
            inf 37471.86283972            inf 22277.897611
            inf  3250.62924497            inf            inf
            inf            inf            inf            inf
            inf            inf            inf            inf
            inf]
[8.48776799e+02            inf 5.87005090e+03            inf
 3.03793615e+04            inf 1.12910734e+05            inf
 2.32123477e+05            inf 2.39378418e+05            inf
 8.47638456e+04            inf 1.91178743e+03 4.84127450e+01
 4.84682588e+01 4.85294076e+01 4.85958162e+01 4.86671141e+01
 4.87429484e+01 4.88229917e+01 4.89069402e+01 4.89945167e+01
 4.90854675e+01]

It gives the attached plot.

[jackaraz]

Hi @APMDSLHC, that is because I set it to inf if it can not converge to a proper value within a given number of iterations. Could you try increasing the iteration_threshold? If that doesn't fix it we need to think of a better solution.

[jackaraz]

Actually, I think I found the problem Tim; please try with the latest push.

[APMDSLHC]

Yep, seems to work fine, thanks.

[jackaraz]

Thanks Tim, closing the issue.