Sample from the joint distribution

In the previous section (Combine independent “analyses” into a joint distribution), having fitted our joint PDF to a single sample, we obtained a second JointDistribution object as output, whose parameters are set to their maximum likelihood estimators given that sample. Its probability distribution is therefore now fully-specified, and we can sample from it:

>>> samples = joint_fitted.sample(3)
>>> print(to_numpy(samples))
{'Test normal::x_theta': array([0., 0., 0.], dtype=float32),
 'Test normal::x': array([10.776527 , 10.519511 ,  9.5911875], dtype=float32),
 'Test binned::x': array([[[-3.0337536, -1.6310872]],
                          [[-4.0343146,  1.5532861]],
                          [[-2.800459 , -0.7948484]]], dtype=float32),
 'Test binned::n': array([[[ 8., 40.]],
                          [[ 2., 50.]],
                          [[ 7., 60.]]], dtype=float32)}

As before, we can fit our distribution to these samples, and this time we will obtain maximum likelihood estimators for each sample independently:

>>> q_3, joint_fitted_3, par_dict_3 = joint.fit_all(samples)
>>> print(to_numpy(par_dict_3["all"]))
{'Test normal':
 {'mu': array([10.776527 , 10.519511 ,  9.5911875], dtype=float32),
  'theta': array([0., 0., 0.], dtype=float32),
  'sigma_t': 0.0},
 'Test binned':
 {'s': array([[[ 0.5640617 , -1.586598  ]],
              [[-0.82253313, -0.7777322 ]],
              [[ 0.22200736,  0.68772215]]], dtype=float32),
  'theta': array([[[-1.5168768 , -0.4077718 ]],
                  [[-2.0171573 ,  0.38832152]],
                  [[-1.4002295 , -0.1987121 ]]], dtype=float32)}}

Note that sigma_t is not considered a free parameter in the BinnedAnalysis class, which is why it still only has one value (for more on this see the BinnedAnalysis class documentation).

Also note that we cannot fit to all the samples simultaneously, i.e. we don’t use each sample as independent information all contributing simultaneously to knowledge of the underlying parameters. JMCTF is designed for performing Monte Carlo simulations of scientific experiments that run just once (such as a search for new particles at the Large Hadron Collider), so each sample is treated as pseudodata whose main purpose is to help us understand the distributions of test statistics. In this view each sample is an independent pseudo-experiment. If an experiment is in fact to be run multiple times in reality, then the PDF of the underlying Analysis class needs to reflect this by using random variables of the appropriate dimension; or for example by using two normal random variables rather than one if the experiment runs twice.

But back to the example. What we did here was a little weird; we sampled from a distribution that was itself fit to some other samples. More usually, we would sample from a distribution with parameters fixed to some “null hypothesis” (or “alternate hypothesis”) values, based on some theoretical motivation. To do this, we can either create the original JointDistribution object with fixed parameters, or fix them in an existing JointDistribution object. But to do this, we need to understand the parameter structure expected by our object. This can be introspected using the get_parameter_structure method:

>>> free, fixed, nuis = joint.get_parameter_structure()
>>> print("free:", free)
>>> print("fixed:", fixed)
>>> print("nuis:", nuis)
free: {'Test normal': {'mu': 1}, 'Test binned': {'s': 2}}
fixed: {'Test normal': {'sigma_t': 1}, 'Test binned': {}}
nuis {'Test normal': {'theta': 1}, 'Test binned': {'theta': 2}}

where ‘free’ are the free parameters in each analysis that can be fit to data, ‘fixed’ are non-fittable parameters that must be chosen by the user (or by theory), and ‘nuis’ are nuisance parameters that enter into the fit but which don’t need to be specified when fixing parameters for sampling (“nominal” values for them can be automatically chosen internally).

To fix some “null hypothesis” input parameters for the analyses we need to provide the parameters from the “free” and “fixed” dictionaries above, but not the “nuis” parameters. For example:

>>> null = {'Test normal': {'mu': [0.], 'sigma_t': [1.]}, 'Test binned': {'s': [(0., 0.)]}}
>>> joint_null = joint.fix_parameters(null)

We can then sample from this null distribution and fit the free parameters:

>>> samples = joint_null.sample(1e6)
>>> q_full, joint_fitted_full, par_dicts_full = joint.fit_all(samples,null)
>>> q_nuis, joint_fitted_nuis, par_dicts_nuis = joint.fit_nuisance(samples, null)

As before (i.e. in the Combine independent “analyses” into a joint distribution section), we use the fit_all method of our JointDistribution object to fit all the free parameters to our samples. Last time we didn’t use the second argument, but it is used to supply a parameter dictionary of the “fixed” parameter values to use during the fitting, if any exist. Here we just gave it the full “null hypothesis” parameter dictionary, but all the parameters aside from the fixed parameter ‘sigma_t’ were ignored.

After this, we used the fit_nuisance method to fit only the nuisance parameters in the analyses to data, with all the other “free” parameters fixed to values specified in the parameter dictionary given in the second argument. So here the “null hypothesis” values of ‘mu’ and ‘s’ were taken as fixed, while both ‘theta’ parameters have been fit.

In the next section, Find MLEs and conduct simple likelihood ratio tests, we will see how the results of these fits can be understood in terms of likelihood ratio tests.