negative_yield3.py

import pyhf
import numpy as np

spec = {
    "channels": [
        {
            "name": "ch",
            "samples": [
                {
                    "data": [100.0,100.0],
                    "modifiers": [
                        {"data": None, "name": "mu_sig", "type": "normfactor"},
                        {"data": {"hi_data": [200,0], "lo_data": [0,200]}, "name": "np_1", "type": "histosys"},
                        {"data": {"hi_data": [200,0], "lo_data": [0,200]}, "name": "np_2", "type": "histosys"},
                        # normsys setup: uncomment the two lines below, comment out the two lines above
                        {"data": {"hi": 1.10, "lo": 0.90}, "name": "np_1", "type": "normsys"},
                        {"data": {"hi": 1.10, "lo": 0.90}, "name": "np_2", "type": "normsys"},
                    ],
                    "name": "signal",
                },
                {
                    "data": [100.0,100.0],
                    "modifiers": [
                        #{"data": None, "name": "mu_sig", "type": "normfactor"},
                        #{"data": {"hi_data": [175], "lo_data": [25.0]}, "name": "np_1", "type": "histosys"},
                        #{"data": {"hi_data": [175], "lo_data": [25.0]}, "name": "np_2", "type": "histosys"},
                        # normsys setup: uncomment the two lines below, comment out the two lines above
                        {"data": {"hi": 1.01, "lo": 0.99}, "name": "np_bkg_1", "type": "normsys"},
                        #{"data": {"hi": 1.75, "lo": 0.25}, "name": "np_2", "type": "normsys"},
                    ],
                    "name": "background",
                }
            ],
        }
    ],
    "measurements": [{"config": {"parameters": [], "poi": "mu_sig"}, "name": "meas"}],
    "observations": [{"data": [115,90], "name": "ch"}],
    "version": "1.0.0",
}


spec2 = {
    "channels": [
        {
            "name": "ch",
            "samples": [
                {
                    "data": [100.0,100.0],
                    "modifiers": [
                        {"data": None, "name": "mu_sig", "type": "normfactor"},
                        {"data": {"hi_data": [220,0], "lo_data": [0,180]}, "name": "np_1", "type": "histosys"},
                        {"data": {"hi_data": [220,0], "lo_data": [0,180]}, "name": "np_2", "type": "histosys"},
                        # normsys setup: uncomment the two lines below, comment out the two lines above
                        #{"data": {"hi": 1.10, "lo": 0.90}, "name": "np_1", "type": "normsys"},
                        #{"data": {"hi": 1.10, "lo": 0.90}, "name": "np_2", "type": "normsys"},
                    ],
                    "name": "signal",
                },
                {
                    "data": [100.0,100.0],
                    "modifiers": [
                        #{"data": None, "name": "mu_sig", "type": "normfactor"},
                        #{"data": {"hi_data": [175], "lo_data": [25.0]}, "name": "np_1", "type": "histosys"},
                        #{"data": {"hi_data": [175], "lo_data": [25.0]}, "name": "np_2", "type": "histosys"},
                        # normsys setup: uncomment the two lines below, comment out the two lines above
                        {"data": {"hi": 1.01, "lo": 0.99}, "name": "np_bkg_1", "type": "normsys"},
                        #{"data": {"hi": 1.75, "lo": 0.25}, "name": "np_2", "type": "normsys"},
                    ],
                    "name": "background",
                }
            ],
        }
    ],
    "measurements": [{"config": {"parameters": [], "poi": "mu_sig"}, "name": "meas"}],
    "observations": [{"data": [115,90], "name": "ch"}],
    "version": "1.0.0",
}

def inspect_it(spec):
    ws = pyhf.Workspace(spec)
    model = ws.model()
    data = ws.data(model)

    print("data ", data)
    print(f"model.config.par_order: {model.config.par_order}")

    bestfit_pars, twice_nll = pyhf.infer.mle.fit(data, model, return_fitted_val=True)

    print(f"bestfit_pars {bestfit_pars}")
    bf_mu0 = np.asarray(bestfit_pars).copy()
    bf_mu0[model.config.poi_index] = 0

    yields = model.expected_actualdata(model.config.suggested_init())
    yields_mle = model.expected_actualdata(bestfit_pars)
    yields_poiZero = model.expected_actualdata(bf_mu0)

    print(f"yields               {yields}",)
    print(f"background yields    {yields_poiZero}")
    print(f"signal yields        {yields_mle}")
    print(f"difference in yields {yields_mle - yields_poiZero}")
    print("-"*20)

inspect_it(spec)
inspect_it(spec2)