1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
# author : S. Mandalia
# s.p.mandalia@qmul.ac.uk
#
# date : April 04, 2018
"""
Likelihood functions for the BSM flavor ratio analysis
"""
from __future__ import absolute_import, division, print_function
from copy import deepcopy
from functools import partial
import numpy as np
import scipy
from scipy.stats import multivariate_normal, truncnorm
from golemflavor import fr as fr_utils
from golemflavor import gf as gf_utils
from golemflavor.enums import Likelihood, ParamTag, PriorsCateg, StatCateg
from golemflavor.misc import enum_parse, gen_identifier, parse_bool
def GaussianBoundedRV(loc=0., sigma=1., lower=-np.inf, upper=np.inf):
"""Normalized Gaussian bounded between lower and upper values"""
low, up = (lower - loc) / sigma, (upper - loc) / sigma
g = scipy.stats.truncnorm(loc=loc, scale=sigma, a=low, b=up)
return g
def multi_gaussian(fr, fr_bf, smearing, offset=-320):
"""
Multivariate Gaussian log likelihood.
Parameters
----------
fr : List[float], length 3
The flavour composition to evaluate at.
fr_bf : List[float], length 3
The bestfit / injected flavour composition.
smearing : float
The amount of smearing.
offset : float, optional
An amount to offset the magnitude of the log likelihood.
Returns
----------
llh : float
The log likelihood evaluated at `fr`.
"""
cov_fr = np.identity(3) * pow(smearing, 2)
return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset
def llh_argparse(parser):
parser.add_argument(
'--stat-method', default='bayesian',
type=partial(enum_parse, c=StatCateg), choices=StatCateg,
help='Statistical method to employ'
)
def lnprior(theta, paramset):
"""Priors on theta."""
if len(theta) != len(paramset):
raise AssertionError(
'Length of MCMC scan is not the same as the input '
'params\ntheta={0}\nparamset={1}'.format(theta, paramset)
)
for idx, param in enumerate(paramset):
param.value = theta[idx]
ranges = paramset.ranges
for value, range in zip(theta, ranges):
if range[0] <= value <= range[1]:
pass
else: return -np.inf
prior = 0
for param in paramset:
if param.prior is PriorsCateg.GAUSSIAN:
prior += GaussianBoundedRV(
loc=param.nominal_value, sigma=param.std
).logpdf(param.value)
elif param.prior is PriorsCateg.LIMITEDGAUSS:
prior += GaussianBoundedRV(
loc=param.nominal_value, sigma=param.std,
lower=param.ranges[0], upper=param.ranges[1]
).logpdf(param.value)
return prior
def triangle_llh(theta, args, asimov_paramset, llh_paramset):
"""Log likelihood function for a given theta."""
if len(theta) != len(llh_paramset):
raise AssertionError(
'Length of MCMC scan is not the same as the input '
'params\ntheta={0}\nparamset]{1}'.format(theta, llh_paramset)
)
hypo_paramset = asimov_paramset
for param in llh_paramset.from_tag(ParamTag.NUISANCE):
hypo_paramset[param.name].value = param.value
spectral_index = -hypo_paramset['astroDeltaGamma'].value
# Assigning llh_paramset values from theta happens in this function.
fr = fr_utils.flux_averaged_BSMu(theta, args, spectral_index, llh_paramset)
flavor_angles = fr_utils.fr_to_angles(fr)
# print('flavor_angles', list(map(float, flavor_angles)))
for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
param.value = flavor_angles[idx]
if args.likelihood is Likelihood.GOLEMFIT:
llh = gf_utils.get_llh(hypo_paramset)
elif args.likelihood is Likelihood.GF_FREQ:
llh = gf_utils.get_llh_freq(hypo_paramset)
return llh
def ln_prob(theta, args, asimov_paramset, llh_paramset):
dc_asimov_paramset = deepcopy(asimov_paramset)
dc_llh_paramset = deepcopy(llh_paramset)
lp = lnprior(theta, paramset=dc_llh_paramset)
if not np.isfinite(lp):
return -np.inf
return lp + triangle_llh(
theta, args=args, asimov_paramset=dc_asimov_paramset,
llh_paramset=dc_llh_paramset
)
|
