pyhf and TRexFitter errors missmatching

[efemberg11]

Hello everyone!
I'm learning how to work with pyhf and trying to compare a fit results (bestfit and errors) of pyhf (0.7.0) and trexf (TRexFitter), using one workspace produced by trexf.

fit values are consistent for both softs, but errors are missmatched.

mu_XS_VVV   1.0007162260432854 (bestfit)   1.139219336823854 (error)             (pyhf)
mu_XS_VVV          1 (bestfit)      +/-  1.35517   (limited)  (error)      (trexf)

Fit provided for the same nuisance params, and same missmatching for data fit case

mu_XS_VVV   -0.251121209117531(bestfit)   0.9788192449023873 (error)
mu_XS_VVV	  = -0.248594	(bestfit)                +/-  1.20699 (error)	(limited)

(another channel)

mu_XS_VVV   2.414723401131362(bestfit)   1.4392019056316985 (error)
mu_XS_VVV         = 2.46056  (bestfit)    +/-  3.33896 (error)  (limited)

I have no idea what to do with it, or is it normal?
Please, can anyone advise what to do and where to go?
With regards, Ramdas.

also my pyhf code

def make_model(channel_list):
    spec["channels"] = [c for c in spec["channels"] if c["name"] in channel_list]

    w = pyhf.Workspace(spec)
    wstat = w
    m = wstat.model(
        measurement_name="OBS_FULL_FIT_VVV",
        modifier_settings={
            "normsys": {"interpcode": "code4"},
            "histosys": {"interpcode": "code4p"},

        },
    )
    d = wstat.data(m)
    return wstat, m, d

def fitresults(chanel, constraints=None):
    _, model, data = make_model(chanel)

    pyhf.set_backend("numpy", pyhf.optimize.minuit_optimizer(verbose=True, errordef=0.5, tolerance=0.01))

    constraints = constraints or []
    init_pars = model.config.suggested_init()
    fixed_params = model.config.suggested_fixed()
    bounds = model.config.suggested_bounds()
    bounds[model.config.poi_index] = [-5, 5]
    for idx, fixed_val in constraints:
        init_pars[idx] = fixed_val
        fixed_params[idx] = True

    result = pyhf.infer.mle.fit(
        data,
        model,
        init_pars=init_pars,
        fixed_params=fixed_params,
        par_bounds=bounds,
        return_uncertainties=True,
    )
    bestfit = result[:, 0]
    errors = result[:, 1]
    return model, data, bestfit, errors

def asimov_fitresults(chanel, constraints=None):
    _, model, data = make_model(chanel)

    pyhf.set_backend("numpy", pyhf.optimize.minuit_optimizer(verbose=True, errordef=0.5, tolerance=0.01))

    constraints = constraints or []
    init_pars = model.config.suggested_init()
    fixed_params = model.config.suggested_fixed()
    bounds = model.config.suggested_bounds()
    bounds[model.config.poi_index] = [-5, 5]
    for idx, fixed_val in constraints:
        init_pars[idx] = fixed_val
        fixed_params[idx] = True

    mu_test = 1
    as_data = pyhf.infer.calculators.generate_asimov_data(
        mu_test,
        data,
        model,
        init_pars=init_pars,
        fixed_params=fixed_params,
        par_bounds=bounds)

    result = pyhf.infer.mle.fit(
        as_data,
        model,
        init_pars=init_pars,
        fixed_params=fixed_params,
        par_bounds=bounds,
        return_uncertainties=True,
    )
    bestfit = result[:, 0]
    errors = result[:, 1]
    return model, data, bestfit, errors

spec = pyhf.readxml.parse('OBS_FULL_FIT_VVV/RooStats/OBS_FULL_FIT_VVV.xml', Path.cwd())
model,data,bestfit,errors = fitresults('three_lep_presel_atLeast_3jets')
print(model.config.channels, " ", model.config.poi_name, " ", bestfit[model.config.poi_index], " ", errors[model.config.poi_index])

[matthewfeickert]

fit values are consistent for both [fits], but errors are missmatched.

@efemberg11 I don't think we can say much without having access to the workspaces and to your TRexFitter code. Have you tried using different values of tolerance?

If you were able to show agreement in an older version of pyhf and iminuit but disagreement now with pyhf v0.7.0 and iminuit v2.17.0 then the pyhf dev team can try to investigate, but without more information I'm not sure what you would like for us to try to do.

[efemberg11]

I tried to change tolerance, there was very small changing at results.
Here I made repository with my codes and workspace.
https://github.com/efemberg11/pyhf-Trexfitter_problem

[alexander-held]

Hi @Efemberg, there are a few things here that contribute to the difference you see:

• The errordef parameter should be set to 1.0, as the function being minimized is currently twice the NLL (change to standard NLL fit objective #1006).
• You are currently using MINOS for the POI uncertainty in TRExFitter but not pyhf.
• The POI bounds you have are [-7, 7] for TRExFitter and [-5, 5] for pyhf (not far enough away from the best-fit points, so this matters a bit here).

With a consistent errordef, bounds of [-7, 7] and without MINOS the results match:

• Asimov
  • TRExFitter: 1 +/- 1.47738
  • pyhf: 1.000000013607906 +/- 1.4773487926962865
• Data
  • TRExFitter: -2.80028 +/- 1.43134
  • pyhf: -2.8002789622557653 +/- 1.4312962850983402

Thanks for sharing the setup in https://github.com/efemberg11/pyhf-Trexfitter_problem! For future reference, the replacement file (my_Common_XS_unc.txt) is missing (but none of the associated uncertainties are used, so they can be commented out). You also seem to have some custom options that you presumably implemented in your version of TRExFitter.

[efemberg11]

Giant thank you! Now I understand. I thought that we use standart NLL.
By the way about MINOS, have you any examples how to calculate minos upper and lower errors in pyhf?

[alexander-held]

You can use cabinetry which supports this via a keyword argument in cabinetry.fit.fit. See the implementation, essentially you'll want to use return_result_obj=True in pyhf.infer.mle.fit and call .minuit.minos() on that.

[alexander-held]

Here is a standalone example: https://gist.github.com/alexander-held/1d45846ed51917d0de895be22b25e826.