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.