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
|
# author : S. Mandalia
# s.p.mandalia@qmul.ac.uk
#
# date : April 04, 2018
"""
Likelihood functions for the BSM flavour ratio analysis
"""
from __future__ import absolute_import, division
import argparse
from functools import partial
import numpy as np
from scipy.stats import multivariate_normal
import GolemFitPy as gf
from utils import fr as fr_utils
from utils import gf as gf_utils
from utils.enums import Likelihood, ParamTag
from utils.misc import enum_parse
def gaussian_llh(fr, fr_bf, sigma):
"""Multivariate gaussian likelihood."""
cov_fr = np.identity(3) * sigma
return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr))
def likelihood_argparse(parser):
parser.add_argument(
'--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood),
choices=Likelihood, help='likelihood contour'
)
def lnprior(theta, paramset):
"""Priors on theta."""
ranges = paramset.ranges
for value, range in zip(theta, ranges):
if range[0] <= value <= range[1]:
pass
else: return -np.inf
return 0.
def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
"""-Log likelihood function for a given theta."""
if len(theta) != len(mcmc_paramset):
raise AssertionError(
'Length of MCMC scan is not the same as the input '
'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, mcmc_paramset)
)
for idx, param in enumerate(mcmc_paramset):
param.value = theta[idx]
hypo_paramset = asimov_paramset
for param in mcmc_paramset.from_tag(ParamTag.NUISANCE):
hypo_paramset[param.name].value = param.value
if args.fix_source_ratio:
fr1, fr2, fr3 = args.source_ratio
else:
fr1, fr2, fr3 = fr_utils.angles_to_fr(
mcmc_paramset.from_tag(ParamTag.SRCANGLES, values=True)
)
bsm_angles = mcmc_paramset.from_tag(
[ParamTag.SCALE, ParamTag.MMANGLES], values=True
)
u = fr_utils.params_to_BSMu(
theta = bsm_angles,
dim = args.dimension,
energy = args.energy,
no_bsm = args.no_bsm,
fix_mixing = args.fix_mixing,
fix_scale = args.fix_scale,
scale = args.scale
)
fr = fr_utils.u_to_fr((fr1, fr2, fr3), u)
for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
param.value = fr[idx]
# print 'hypo_paramset', hypo_paramset
if args.likelihood is Likelihood.FLAT:
return 1.
elif args.likelihood is Likelihood.GAUSSIAN:
fr_bf = args.measured_ratio
return gaussian_llh(fr, fr_bf, args.sigma_ratio)
elif args.likelihood is Likelihood.GOLEMFIT:
return gf_utils.get_llh(fitter, hypo_paramset)
def ln_prob(theta, args, fitter, asimov_paramset, mcmc_paramset):
lp = lnprior(theta, paramset=mcmc_paramset)
if not np.isfinite(lp):
return -np.inf
return lp + triangle_llh(
theta, args=args, asimov_paramset=asimov_paramset,
mcmc_paramset=mcmc_paramset, fitter=fitter
)
|
