StatError¶

In [2]:

import pyhf
import importlib

In [48]:

spec = {
    'channels': [
        {
            'name': 'firstchannel',
            'samples': [
                {
                    'name': 'mu',
                    'data': [10.0, 10.0],
                    'modifiers': [{'name': 'mu', 'type': 'normfactor', 'data': None}],
                },
                {
                    'name': 'bkg1',
                    'data': [50.0, 70.0],
                    'modifiers': [
                        {
                            'name': 'stat_firstchannel',
                            'type': 'staterror',
                            'data': [10.0, 10.0],
                        }
                    ],
                },
                {
                    'name': 'bkg2',
                    'data': [30.0, 20.0],
                    'modifiers': [
                        {
                            'name': 'stat_firstchannel',
                            'type': 'staterror',
                            'data': [5.0, 5.0],
                        }
                    ],
                },
                {'name': 'bkg3', 'data': [20.0, 15.0], 'modifiers': []},
            ],
        },
        #         {
        #             'name': 'secondchannel',
        #             'samples': [
        #                 {
        #                     'name': 'bkg2',
        #                     'data': [30.0],
        #                     'modifiers': [
        #                         {'name': 'stat_secondchannel', 'type': 'staterror', 'data': [5.]}
        #                     ]
        #                 }
        #             ]
        #         }
    ]
}
p = pyhf.Model(spec)

In [49]:

# se = p.config.modifier('stat_firstchannel')

In [50]:

tensorlib = pyhf.tensorlib

In [51]:

inquad = tensorlib.sqrt(tensorlib.sum(tensorlib.power(se.uncertainties, 2), axis=0))
totals = tensorlib.sum(se.nominal_counts, axis=0)
uncrts = tensorlib.divide(inquad, totals)
uncrts

Out[51]:

array([0.16666667, 0.25      ])

In [52]:

se.pdf(se.auxdata, [1.0, 1.0])

Out[52]:

array([2.39365368, 1.59576912])

In [53]:

# p.config.par_slice('stat_firstchannel')

In [54]:

p.config.auxdata

Out[54]:

[1.0, 1.0]

In [55]:

p.spec['channels'][0]['samples'][1]

Out[55]:

{'name': 'bkg1',
 'data': [50.0, 70.0],
 'modifiers': [{'name': 'stat_firstchannel',
   'type': 'staterror',
   'data': [10.0, 10.0]}]}

In [56]:

import pyhf
import json
import logging
from pyhf import runOnePoint, Model
from pyhf.simplemodels import uncorrelated_background


def invert_interval(testmus, cls_obs, cls_exp, test_size=0.05):
    point05cross = {'exp': [], 'obs': None}
    for cls_exp_sigma in cls_exp:
        yvals = [x for x in cls_exp_sigma]
        point05cross['exp'].append(
            np.interp(test_size, list(reversed(yvals)), list(reversed(testmus)))
        )

    yvals = cls_obs
    point05cross['obs'] = np.interp(
        test_size, list(reversed(yvals)), list(reversed(testmus))
    )
    return point05cross


def plot_results(testmus, cls_obs, cls_exp, test_size=0.05):
    plt.plot(mutests, cls_obs, c='k')
    for i, c in zip(range(5), ['grey', 'grey', 'grey', 'grey', 'grey']):
        plt.plot(mutests, cls_exp[i], c=c)
    plt.plot(testmus, [test_size] * len(testmus), c='r')
    plt.ylim(0, 1)

In [57]:

pdf = p
%pylab inline
data = [100.0, 100.0] + pdf.config.auxdata

init_pars = pdf.config.suggested_init()
par_bounds = pdf.config.suggested_bounds()

mutests = np.linspace(0, 5, 41)
tests = [
    runOnePoint(muTest, data, pdf, init_pars, par_bounds)[-2:] for muTest in mutests
]
cls_obs = [test[0] for test in tests]
cls_exp = [[test[1][i] for test in tests] for i in range(5)]

plot_results(mutests, cls_obs, cls_exp)

Populating the interactive namespace from numpy and matplotlib