How to draw toy MC?

[zjkjsd]

Dear experts,

Could you please assist me in understanding how to use pyhf for drawing toy Monte Carlo (MC) samples? I came across the following code snippet, which seems to draw toys from the template model:

toys = model.make_pdf(pyhf.tensorlib.astensor(np.asarray(model.config.suggested_init()))).sample((n_toys,))

1. My understanding is that we first sum all the components and then draw samples from the Poisson distribution defined for each bin, with the mean set to the count in the corresponding template bin. Is this understanding correct? And is this the correct approach for drawing toys?

2. Furthermore, if I wish to draw toys with scaling factors applied to the signal and background components (for generating validation plots like linearity plots), should I replace model.config.suggested_init() with a list of the scaling factors?

3. If my understanding is correct, would the following code also achieve the desired outcome? (I included it to confirm my understanding.)

@joblib.delayed
def compute_toys(sig_hist, bkg_hist, sig_scale, bkg_scale, itry):
    MC = sig_scale * sig_hist + bkg_scale * bkg_hist
    rng = np.random.default_rng(itry)
    toy = rng.poisson(MC)
    return toy
ntoys=20000
toys = joblib.Parallel(-1)(compute_toys(sig_hist, bkg_hist, 0.5, 0.8, i) for i in range(ntoys))

Thank you for your support!
Boyang

[kratsg]

If the question is about how to throw toys using pyhf, we have calculators that do this for you. In particular, see this example in the documentation as well as this example reproducing a key plot in the HistFactory paper.

If you're asking about how to do it yourself, look here:

pyhf/src/pyhf/infer/calculators.py

Lines 756 to 779 in c096b33

	tensorlib, _ = get_backend()
	sample_shape = (self.ntoys,)
	signal_pars = fixed_poi_fit(
	poi_test,
	self.data,
	self.pdf,
	self.init_pars,
	self.par_bounds,
	self.fixed_params,
	)
	signal_pdf = self.pdf.make_pdf(signal_pars)
	signal_sample = signal_pdf.sample(sample_shape)
	bkg_pars = fixed_poi_fit(
	1.0 if self.test_stat == 'q0' else 0.0,
	self.data,
	self.pdf,
	self.init_pars,
	self.par_bounds,
	self.fixed_params,
	)
	bkg_pdf = self.pdf.make_pdf(bkg_pars)
	bkg_sample = bkg_pdf.sample(sample_shape)

. In particular, you need to evaluate the model at the fitted MLE and use the $\hat{\theta}$ parameters to initialize your pdf to sample from (rather than the default suggested_init). At that point, we just sample as you typically expect from the "scipy API" in the same way.

1. Yes, but it's Poisson and Gaussian.
2. Yes. But see my description/explanation above.
3. That code seems to be more generic about using a Poisson to make toys.

[alexander-held]

The key thing missing in your example @zjkjsd are the fluctuations of auxiliary data ("global observables"). For nuisance parameters with Gaussian constraints, the corresponding auxiliary data is drawn from Gaussians (that's where the Gaussian in @kratsg's 1. appears). If you look at the output of pyhf's implementation, you will see that both the counts per bin and the auxiliary data fluctuate.

[zjkjsd]

Hi Giordon and Alex,

Thank you for your help, and I have a question for the code below:

# define samples and parameters
signal = 10000   # single bin
background = 0 # turned off for simplicity
background_uncertainty = 0
model = pyhf.simplemodels.uncorrelated_background(
    [signal], [background], [background_uncertainty]
)
# define toy par and generate toys
toy_pars = [0.1, 0] 
pdf_toy = model.make_pdf(pyhf.tensorlib.astensor(toy_pars))
toys = pdf_toy.sample((1000,))

# calculate the standard deviation of the actual data over all toys
sd = np.std(toys[:,0])
print(sd)

I would expect the standard deviation sd to be around 0.1 * sqrt(10000)=10, but the result from the code above is 31.6, why is that? It seems to be doing sqrt( 0.1 * 10000), but I thought this was incorrect.

[zjkjsd]

Thank you, Alex! Yeah, I don't use observed data in the linearity test. Could you please confirm the steps below are reasonable procedures according to our discussions?

1. run toy_pars = asimov_parameters(model)
2. change the normfactor parameters according to my linearity test, e.g. toy_pars[:2] = [0.1, 0.1]
3. generate toys with model.make_pdf(pyhf.tensorlib.astensor(toy_pars)).sample((ntoys,))

I am trying to understand why asimov_parameters should be used for my linearity test, i.e. why we want to set other parameters to values that maximize the associated constraint terms. Now I guess that: varying the signal strength does not change our understanding of the systematics, so we should use our best knowledge in the constraint terms, which correspond to 1 for shapesys and 0 for normsys/histosys. Is this correct? Thanks again!

[alexander-held]

I do not know what you are ultimately trying to do and why you are generating toys. I expect you have a specific reason to generate them and in that case the choice of parameters would be motivated by what you are aiming to do with the toys. Generically all I can say is that in the absence of any additional knowledge (like from incorporating information from observed data), the most sensible setting for the parameters is the one maximizing the constraint terms and that holds true independent of any normalization factor values as you say.

I am not sure what linearity you are looking for and why you need toys for that. If you want to make sure you fit the injected normalization factor value again then this is something that happens by design (in the absence of technical problems) and it is unclear to me which additional information toys provide.

[zjkjsd]

Thank you for the explanation! Again, this is very helpful! I have 2 goals for generating toys:

1. fit the injected normalization factor value
2. calculate the uncertainty of the fitted normalization factors

It's good to know some happen by design, thanks! I am trying to do this: I have 2 workspaces whose background samples are constructed from different regions, w_sig from the signal region and w_sb from sidebands. I am trying to generate toys with my w_sig workspace and fit toys with w_sb. This way, I can estimate how well the sidebands describe the signal region, and if the sidebands are good enough (i.e. doesn't bias the signal strength). Does it make sense to use toys for my goals?

[alexander-held]

In general I would say that 1) happens by design and for 2) you can also look at the uncertainty in your fit to the Asimov dataset which might already tell you what you are after. See https://arxiv.org/abs/1007.1727 for more about this. We commonly use toys in cases where the asymptotic results derived in that paper do not hold. Another common use case is the estimation of trial factors for the look elsewhere effect. There are applications beyond that but I would say they are more specific and usually more rare.

I probably do not understand the setup you have. What sometimes happens is that the background estimate for a signal region is obtained from sidebands which you might fit first in a separate step or simultaneously with your signal region. I do not spontaneously see a need for toys in such a setup. I would recommend chatting with colleagues from Belle II which might be more familiar with your setup and can perhaps suggest common approaches used for such cases.

[zjkjsd]

Thanks for the explanation and suggestions! I will talk with my advisor and other Belle II people. To clarify my use of sidebands, I am taking the sidebands as a single sample to model a bkg component (FakeD) in the signal region, this way I don't need to use MC for this bkg (potential Data/MC discrepancy). I have heard fitting signal region and sidebands simultaneously, but I am not sure how to do it for my analysis. It sounds like I need to have templates for FakeD in both regions, and I am not sure how to avoid using MC.