David Rousso toys notebook converted to script

roussou_debug.py

import numpy as np
import pyhf

stopRHadron_spec = {
    'channels': [
        {
            'name': 'channel1',
            'samples': [
                {
                    'data': [1.364790054231882],
                    'modifiers': [
                        {'data': None, 'name': 'lumi', 'type': 'lumi'},
                        {'data': None, 'name': 'mu_Sig', 'type': 'normfactor'},
                        {
                            'data': {'hi': 1.228925751097454, 'lo': 0.7710742489025461},
                            'name': 'ucs_StopRHadron_1100_100000',
                            'type': 'normsys',
                        },
                    ],
                    'name': 'StopRHadron_1100_100000',
                },
                {
                    'data': [0.43],
                    'modifiers': [
                        {
                            'data': [0.16],
                            'name': 'staterror_channel1',
                            'type': 'staterror',
                        }
                    ],
                    'name': 'Bkg',
                },
            ],
        }
    ],
    'measurements': [
        {
            'config': {
                'parameters': [
                    {
                        'auxdata': [1.0],
                        'bounds': [[0.5, 1.5]],
                        'inits': [1.0],
                        'name': 'lumi',
                        'sigmas': [0.017],
                    }
                ],
                'poi': 'mu_Sig',
            },
            'name': 'meas',
        }
    ],
    'observations': [{'data': [0.0], 'name': 'channel1'}],
    'version': '1.0.0',
}


def main():
    pyhf.set_backend("jax")

    workspace = pyhf.Workspace(stopRHadron_spec)
    model = workspace.model()
    observations = [0]
    data = observations + model.config.auxdata
    print(f"data: {data}")
    test_poi = 1.0

    print("# Asymptotics\n")
    cls_obs, cls_exp = pyhf.infer.hypotest(
        test_poi, data, model, return_expected_set=True, calctype="asymptotics"
    )
    print(f"cls_obs: {cls_obs}")
    print(f"cls_exp: {cls_exp}")

    n_toys = 10000
    print("\n# Toys")
    cls_obs, cls_exp = pyhf.infer.hypotest(
        test_poi,
        data,
        model,
        return_expected_set=True,
        calctype="toybased",
        ntoys=n_toys,
    )
    print(f"\ncls_obs: {cls_obs}")
    print(f"cls_exp: {cls_exp}")


if __name__ == '__main__':
    main()

Running this produces the following output

data: [0, 1.0, 1.0, 0.0]
# Asymptotics

cls_obs: [0.15568652]
cls_exp: [[0.02969   ]
 [0.08653433]
 [0.22822617]
 [0.49778549]
 [0.8050249 ]]

# Toys
                                                                                                                                                                                                            
cls_obs: [0.24197901]
cls_exp: [[0.25022422]
 [0.26773903]
 [1.        ]
 [1.        ]
 [1.        ]]

Note: very small numbers so asymptotics may not set in yet.
The systematic variations ucs_StopRHadron_1100_100000 are both less than the nominal, which may cause some strange behavior.

Is there a nice way to plot a histogram of the test statistic?

Is there a nice way to plot the values of the parameters being used to generate the toys (vs. mu)?

below I show brazil bands for scans of µ \in [0.1, 5, 50]

asymptotics

toys (1k toys per µ)