Fitting problem: Finding the signal that best fits the data from a set of possible signals

[IzaV10]

Hi all,

I am struggling with implementing the following maximum likelihood fit using pyhf and am wondering if the option exists. I want to find the maximum likelihood for the following setup:

• Let’s say that we have some expected data and a background model, B. There is a set of possible signals, S(C), which depend on C.

• The signal, S(C), for any C, can be obtained as S(C)=S(C1)*coeff_1(C)+ S(C2)*coeff_2(C)+ S(C3)*coeff_3(C), where S(C_x) is the expected signal for C_x and coeff_x(C) are coefficients.

• The goal is to find which S(C)+B model fits best the expected data. In other words, which value of C maximizes the likelihood.

There are two ways (that I can think of) how to do this:

1. The POI is an equation that depends on C and a few templets of S(C_x) [POI(C)=S(C1)*coeff_1(C)+ S(C2)*coeff_2(C)+ S(C3)*coeff_3(C)].

2. Use S(C)=S(C1)*coeff_1(C)+ S(C2)*coeff_2(C) to generate signal histograms for several values of C, so the POI=C and find by brute force which signal fits the expected data best.

Is any of the above mentioned methods supported by pyhf?
Thank you very much for your time and help.

Cheers,
Iza

[alexander-held]

Hi Iza, are S(C1), S(C2) and S(C3) three different constant templates, or are C1/C2/C3 also functions of C?

The second approach sounds correct to me. The way to solve this in ROOT HistFactory is via expressions (AddPreprocessFunction). You have three signal templates, and would scale each of them by a function that depends on your POI C. In ROOT there is no need to manually generate many histograms for different C, the three histograms you already have should be enough in principle. This is similar to how you can specify up/down variations for systematic uncertainties, and then the HistFactory interpolation/extrapolation prescription gives you values for any parameter value. In this case the interpolation/extrapolation would be explicitly specified via your formula.

I do not understand the first approach, what you call POI here sounds like the signal yield per bin to me, which is not a single parameter unless you have only a single bin (and if you have multiple bins, the yields are not all independent).

Expressions for normalization factors are not yet supported in pyhf (see #850 and #1188), so there is no easy way to implement this currently. Working around that limitation may be possible with some limitations, although tedious: you could generate workspaces parameterized as a function of C (only the signal template changes across workspaces, but C is not a parameter in any of the workspaces). I assume this is what you mean by "brute force"? Then you could fit each workspace, and picking the best-fit model will tell you what the best-fit C is (the fit of the workspace in which the signal template was built with a specific value of C) and a manual likelihood scan over workspaces will also reveal the uncertainty in C. The problem with that approach is that C is never a parameter in your fit, and you cannot account for possible correlations between C and other parameters in the fit model. That will generally impact the uncertainty calculation. It may instead be easier to give the open PR #1188 a try to implement the correct treatment.

[IzaV10]

Hi,

S(C1), S(C2) and S(C3) are three different constant templates from which any S(C) template can be constructed (you could even construct S(C1) from a linear combination of three other fixed S(C_x) templates for example).

Yes, you suggested "work around" is what I am currently doing.

I had a look at the open PR #1188 but am a bit confused on how to use custom modifier locally? Does on need to copy the whole repository and checkout the branch on which the work has been done?
Can the custom modifier depend on three signals, ie S(C)=S(C1)*coeff_1(C)+ S(C2)*coeff_2(C)+ S(C3)*coeff_3(C)? How would you implement that? Is there any open issues that would need fixing first?

Thanks for your help!

Cheers,
Iza

[alexander-held]

You can install from the branch directly with pip:

pip install git+https://github.com/scikit-hep/pyhf.git@custom_modifiers

Alternatively, you can also clone the repository, check out the branch, and install from within the repository (via pip install .).

@lukasheinrich will know the details better than me, but I believe in the current implementation the modifier is general and you can implement any function you want there, meaning that you could also interpolate between three signals manually.

The related issue #850 deals with a more restricted case, namely writing formulas that depend on normalization factors (like expressions in TRExFitter / AddPreprocessFunction() in ROOT HistFactory). An implementation resolving #850 should be sufficient for your setup.

See also #1188 (comment) and Lukas' reply.

[IzaV10]

Thanks you very much! I will have a go

[IzaV10]

Hi,

When doing so and following the example in #1188 I get the problem described here. @lukasheinrich is the issue mentioned above solvable as well as there maybe some documentation on how to use custom modifier?

Cheers,
Iza