Create an instance of PyFBU
import fbu
myfbu = fbu.PyFBU()
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions. /home/francesco/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`. from ._conv import register_converters as _register_converters
Supply the input distribution to be unfolded as a 1-dimensional list for N bins, with each entry corresponding to the bin content.
myfbu.data = [100,150]
Supply the response matrix where each row corresponds to a truth level bin. The normalization of each row must be the acceptance efficiency of the corresponding bin (e.g. the normalization is 1 for resolution only unfolding). N.B. For now, only squared response matrices are allowed.
myfbu.response = [[0.08,0.02], #first truth bin
[0.02,0.08]] #second truth bin
Define the boundaries of the hyperbox to be sampled for each bin.
myfbu.lower = [0,0]
myfbu.upper = [3000,3000]
Run the MCMC sampling (this step might take up to several minutes for a large number of bins).
myfbu.run()
Assigned Metropolis to truth0 Assigned Metropolis to truth1 100%|██████████| 11000/11000 [00:02<00:00, 3948.16it/s]
Retrieve the N-dimensional posterior distribution in the form of a list of N arrays.
trace = myfbu.trace
print( trace )
[array([915, 814, 775, ..., 975, 850, 850]), array([1523, 1523, 1505, ..., 1770, 1770, 1893])]
Each array corresponds to the projection of the posterior distribution for a given bin.
%matplotlib inline
from matplotlib import pyplot as plt
plt.hist(trace[1],
bins=20,alpha=0.85,
normed=True)
plt.ylabel('probability')
Text(0,0.5,'probability')
One or more backgrounds, with the corresponding normalization uncertainties (gaussian prior), can be taken into account in the unfolding procedure.
myfbu.background = {'bckg1':[20,30],'bckg2':[10,10]}
myfbu.backgroundsyst = {'bckg1':0.5,'bckg2':0.04} #50% normalization uncertainty for bckg1 and 4% normalization uncertainty for bckg2
The background normalization is sampled from a gaussian with the given uncertainty. To fix the background normalization the uncertainty should be set to 0.
Systematic uncertainties affecting signal and background can be taken into account as well with their per-bin relative magnitudes. The prior is gaussian. Each systematics needs to be provided for each background listed at the previous step.
myfbu.objsyst = {
'signal':{'syst1':[0.,0.03],'syst2':[0.,0.01]},
'background':{
'syst1':{'bckg1':[0.,0.],'bckg2':[0.1,0.1]},
'syst2':{'bckg1':[0.,0.01],'bckg2':[0.,0.]}
}
}
Each systematics is treated as fully correlated across signal and the various backgrounds.
The posterior probability for the nuisance parameters is stored in a dictionary of arrays. The correlation among nuisance parameters and with the estimates for the unfolded distribution is preserved in the array ordering.
myfbu.run() #rerun sampling with backgrounds and systematics
Assigned Metropolis to truth0 Assigned Metropolis to truth1 Assigned NUTS to gaus_bckg1_lowerbound__ Assigned NUTS to gaus_bckg2_lowerbound__ Assigned NUTS to gaus_syst1 Assigned NUTS to gaus_syst2 /home/francesco/anaconda3/lib/python3.6/site-packages/pymc3/model.py:384: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`. if not np.issubdtype(var.dtype, float): 100%|█████████▉| 10974/11000 [00:28<00:00, 385.37it/s]/home/francesco/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:467: UserWarning: Chain 0 contains 61 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize. % (self._chain_id, n_diverging)) 100%|██████████| 11000/11000 [00:28<00:00, 385.29it/s]
unfolded_bin1 = myfbu.trace[1]
bckg1 = myfbu.nuisancestrace['bckg1']
plt.hexbin(bckg1,unfolded_bin1,cmap=plt.cm.YlOrRd)
<matplotlib.collections.PolyCollection at 0x7f2616df5d30>