questions about how to use pyhf

[zjkjsd]

Aloha pyhf community,

Thank you all for making pyhf available! My name is Boyang, a physics student at the University of Hawaii. I want to use this package in my inclusive tagging R(D) measurement at Belle II. It would be greatly appreciated if you can help me with the 3 questions below.

1. Does the order of templates in a channel matter?
   I have 10 components in a channel thus 10 templates ( 6 are signal-like, B -> D(*, **) tau(ell) nu; 4 are backgrounds ) in a channel, and only normalizations are fitted. I am wondering if the fitting results would change when varying the order of templates placed in the workspace json file. For example, is it supposed to make any difference if I switch the 1st template to the last one? I guess it would because the model seems to imply that the components are fitted one by one instead of a single multi-dimensional fit (if I understand correctly).

2. is it possible to constrain the ratio of 2 normalizations in the fitting?
   A useful constraint is to set the ratio (B -> D tau nu / B -> D* tau nu) equal to ( B -> D ell nu / B -> D* ell nu) because the 2 components in the 2nd ratio are better fitted. Could you share some insights about how to implement this in pyhf? The fixed_params in pyhf.infer.mle.fit might be relevant, but I am not sure how to use it.

3. Is there a default pyhf way to remove empty bins in 2d fitting?
   I am considering 2d fitting with missing mass squared and momentum of (D+ell, scalar sum). These 2 variables are correlated, so the 2d histogram looks diagonal and many off-diagonal bins are empty. Is there a default way to exclude/remove empty bins in the fitting so that only relevant bins are fitted? The main concern is that the empty bins' padding may confuse the fitter.

I would be really grateful if you can give some feedback on my code below!
Boyang

xedges = np.linspace(-2, 10, 48) # -7.5 for weMiss2, -2 for weMiss3, -2.5 for weMiss4
yedges = np.linspace(0.4, 4.6, 42)
variable_x = 'B0_CMS3_weMissM2'
variable_y = 'p_D_l'

(counts, xedges, yedges) = np.histogram2d(df_merged.query(cut)[variable_x], 
                                              df_merged.query(cut)[variable_y],
                                              bins=[xedges, yedges])
counts = counts.T
    
def fit():
    global parameters, correlations
    parameters, correlations= pyhf.infer.mle.fit(np.concatenate([counts.flat]), model_1, return_uncertainties=True, return_correlations=True)
    sum_event=0
    # calculate the fitted counts and errors
    for i in range(len(parameters)):
        norm = parameters[i,0]
        template_counts = np.sum(model_1.spec['channels'][0]['samples'][i]['data'])
        component_name = model_1.spec['channels'][0]['samples'][i]['name']
        fitted_counts = round(parameters[i,0] * template_counts)
        fitted_error = round(parameters[i,1] * template_counts)
        print(f'\n{component_name} counts:             {fitted_counts}')
        print(f'{component_name} counts uncertainty: {fitted_error}')
        sum_event+= fitted_counts

    print(f'\nTotal fitted counts: {sum_event}')
    print(f'Total fitted counts uncertainty: ??')
    print(f'Total generated counts: {len(df_merged)}')
    print(f'\nFitting parameters: \n{parameters}')
    print(f'Fitting parameter correlations: \n{correlations}')
    print(f'uncertainty: {parameters[0][1]}')

[alexander-held]

Hi @zjkjsd, I'll try to get this started.

1. No, the order does not matter. The caveat to that is that a different template order could in principle lead to a very slightly different result due to floating point math (should be within minimization tolerance). As far as I am aware, pyhf internally orders the templates alphabetical by name so the order in the specification is irrelevant and only the name matters.

the model seems to imply that the components are fitted one by one instead of a single multi-dimensional fit

The fit will involve all the templates you have specified. You could control templates via a normfactor and hold specific templates constant (and set to zero) in the fit, but that is something you would have to configure in your fits and the model itself will by default involve all templates.

2. In general this is addressed by Application of POIs using formulas #850. For the specific case of a ratio, you can actually get this to work without the need for expressions. Here is how you would set this up:

• define normfactor A := B -> D tau nu / B -> D* tau nu,
• define normfactor B := B -> D* tau nu,
• define normfactor C := B -> D* ell nu.
  Now apply the normfactor modifiers:
• sample B -> D tau nu is scaled by both modifiers A and B,
• sample B -> D* tau nu is scaled by B,
• sample B -> D ell nu is scaled by both modifiers A and C (since you set the ratios to be equal),
• sample B -> D* ell nu is scaled by C.

fixed_params is for holding parameters constant in a fit and unrelated to what you are asking for, assuming I understand the question correctly.

3. No, pyhf does not provide any utility for that. (Caveat: there are some things happening internally where pyhf will remove terms in the likelihood corresponding to empty bins, for example for staterror modifiers, perhaps for other cases too). In general it will take all the bins you provide). The HistFactory model employed in pyhf also only handles one-dimensional distributions. You can "unroll" 2d distributions though, and the only place where this can make a difference are modifiers that act across channels (like histosys), which introduce a notion of bins being connected per channel.

Regarding your code snippet:

fitted_error = round(parameters[i,1] * template_counts)

This only is true for the specific case of a single normfactor acting on (a) sample(s). As soon as you have more parameters that can modify the count, you need to consider all of them and their correlations. Instead of extracting numbers from the specification via

np.sum(model_1.spec['channels'][0]['samples'][i]['data'])

you can use the model.expected_data API (or the model.main_model.expected_data version).

print(f'Total fitted counts uncertainty: ??')

For this you need to consider the effects of all parameters (and their correlations) acting on the full model.

I would recommend having a look at cabinetry, which provides these calculations. cabinetry.model_utils.prediction gives you all the calculations at once, while cabinetry.model_utils.yield_stdev does only the uncertainty calculation.

[zjkjsd]

Aloha Alex,

Thank you so much for your response! I understand much more about pyhf now, and your feedback about my code snippet looks very useful! Can I ask some follow-up questions:

2. You understood my question very well! fixed_params should not be needed. I wrote a model below according to your suggestion, does it look ok? Actually, I am not sure what you meant by define and apply normfactor.

model = pyhf.Model(
    {
        'channels': [
            {
                'name': ‘e_lchannel',
                'samples': [
                    {'name': ‘D_tau_nu’, 'data': [10] * 20,
                     'modifiers': [
                            {'name': ‘A’,’type': 'normfactor','data': null},
			    {'name': ‘B’,’type': 'normfactor','data': null},
                        ],
                    },
                    {'name': ‘Dst_tau_nu’, 'data': [10] * 20,
                     'modifiers': [
			    {'name': ‘B’,’type': 'normfactor','data': null},
                        ],
                    },
		   {'name': ‘D_ell_nu’, 'data': [10] * 20,
                     'modifiers': [
                            {'name': ‘A’,’type': 'normfactor','data': null},
			    {'name': ‘C’,’type': 'normfactor','data': null},
                        ],
                    },
		   {'name': ‘Dst_ell_nu’, 'data': [10] * 20,
                     'modifiers': [
			    {'name': ‘C’,’type': 'normfactor','data': null},
                        ],
                    },
                ],
            }
        ]
    },
    poi_name=None, # I was told to set this to None, does it make sense?
)

Does poi_name=None make sense? I am not doing a hypothesis test so I was told to set it as is.
However, I now realized that what I want is really to set the ratio of fitted counts (B -> D tau nu / B -> D* tau nu) equal to ( B -> D ell nu / B -> D* ell nu), not the ratio of normalization. So do I have to use expressions now? If so, should I follow the page and use the function add_custom_modifier? i.e. (Here N stands for normalization factor)

$$ \frac{N(D\tau\nu)\times Template\ counts(D\tau\nu)}{N(D^\ast\tau\nu)\times Template\ counts(D^\ast\tau\nu)} = \frac{N(D\ell\nu)\times Template\ counts(D\ell\nu)}{N(D^\ast\ell\nu)\times Template\ counts(D^\ast\ell\nu)}$$

or

$$ \frac{N(D\tau\nu)}{N(D^\ast\tau\nu)} = R\times\frac{N(D\ell\nu)}{N(D^\ast\ell\nu)}$$

$$ R = \frac{Template\ counts(D^\ast\tau\nu)\times Template\ counts(D\ell\nu)}{Template\ counts(D\tau\nu)\times Template\ counts(D^\ast\ell\nu)}$$

new_params = {
    ‘A’: {'inits': (1.0,), 'bounds': ((0,1),)}, # D*_ell_nu / D_ell_nu
    ‘B’: {'inits': (1.0,), 'bounds': ((0,1),)}, # D_tau_nu
    ‘C’: {'inits': (1.0,), 'bounds': ((0,1),)}, # D_ell_nu
}
R = np.sum(template_1 * template_4 / (template_2 * template_3) )

expanded_pyhf = add_custom_modifier('customfunc', [‘A’,’B’,’C’], new_params)
model = pyhf.Model(
    {
        'channels': [
            {
                'name': ‘e_lchannel',
                'samples': [
                    {'name': ‘Dst_tau_nu’, 'data': [10] * 20,
                     'modifiers': [
                            {'name': ‘D’,'type': 'customfunc','data': {'expr': 'R*A'}},
			    {'name': ‘B’,’type': 'normfactor','data': null},
                        ],
                    },
                    {'name': ‘D_tau_nu’, 'data': [10] * 20,
                     'modifiers': [
			    {'name': ‘B’,’type': 'normfactor','data': null},
                        ],
                    },
		   {'name': ‘Dst_ell_nu’, 'data': [10] * 20,
                     'modifiers': [
                            {'name': ‘A’,’type': 'normfactor','data': null},
			    {'name': ‘C’,’type': 'normfactor','data': null},
                        ],
                    },
		   {'name': ‘D_ell_nu’, 'data': [10] * 20,
                     'modifiers': [
			     {'name': ‘C’,’type': 'normfactor','data': null},
                        ],
                    },
                ],
            }
        ]
    },
    modifier_set=expanded_pyhf,
    poi_name=None,
    validate=False, #??
    batch_size = 1 #??
)

How does this look to you?

[zjkjsd]

3. Do you think padding the empty bins in templates would affect the fitting performance?
   My 2 fitting variables are correlated and more than 500 bins are empty (I have 2000 bins in total) in my 2d histogram. So far, I have been keeping those empty bins in my templates (after flattening the 2d histogram) and as we know they also have to be padded with a small number (otherwise throw an optimization error). Do you think the padding on empty bins will confuse the fitter? If so, I can remove the empty bins before putting them into the templates. If not, I will keep the empty bins, is there a recommended value for padding?

[kratsg]

1/2 look fine. If you're going for a somewhat more complicated expression like that, you will want to use the custom mods approach here.

3. padding won't necessarily affect performance. The idea in pyhf is that the specific bins are masked from the calculation meaning that while they still get included when doing the actual multiplication across, they're dropped/masked at the end when computing the result.