gambit is hosted by Hepforge, IPPP Durham
 GAMBIT  v1.5.0-2191-ga4742ac a Global And Modular Bsm Inference Tool
convolve_with_theory.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 """
3 Convolve chi-squared in a data file with a fractional theory error
4 ==================================================================
5
6 python convolve_with_theory.py <file> <frac-error> <min> <max>
7
8 prints (parameter, convolved chi-squared) from a file containing
9 (parameter, chi-squared).
10 """
11
12 import sys
13 import numpy as np
14 from scipy.interpolate import interp1d
15 from scipy.stats import norm, lognorm
16 from scipy.integrate import quad
17
18
19 def convolve(file_name, frac_error=0.1, min_=0., max_=1., log_normal=True):
20  """
21  Convolve chi-squared in a data file with a fractional theory error
22
23  Args:
24  file_name (str): Data file with columns (parameter, chi-squared).
25  frac_error (float, optional): Fractional theory error on the parameter.
26  min_ (float, optional): Minimum value of parameter.
27  max_ (float, optional): Maximum value of parameter.
28  log_normal (bool, optional): Whether to use log-normal or normal error.
29
30  Returns:
31  list(tuples): List of (parameter, convolved chi-squared)
32  """
33  # Unpack data
34  param, chi_squared = np.loadtxt(file_name, unpack=True)
35
36  # Interpolate likelihood function
37  like = interp1d(param, np.exp(-0.5 * chi_squared),
38  kind='linear', bounds_error=False, fill_value="extrapolate")
39
40  # Make prior for true prediction
41  def prior(x, mu):
42  if log_normal:
43  sigma = np.log(1. + frac_error)
44  dist = lognorm(sigma, scale=mu)
45  else:
46  sigma = frac_error * mu
47  dist = norm(mu, sigma)
48  return dist.pdf(x)
49
50  # Make functions for convolution
51  integrand = lambda x, mu: like(x) * prior(x, mu)
52  convolved = lambda mu: -2. * np.log(
53  quad(integrand, min_, max_, args=(mu))[0]
54  / quad(prior, min_, max_, args=(mu))[0])
55
56  # Perform convolution
57  return [(p, convolved(p)) for p in param]
58
59 if __name__ == "__main__":
60
61  try:
62  FILE_NAME = sys.argv[1]
63  except IndexError:
64  raise RuntimeError(__doc__)
65
66  ARGS = map(float, sys.argv[2:])
67
68  for p, c in convolve(FILE_NAME, *ARGS):
69  print p, c
def convolve(file_name, frac_error=0.1, min_=0., max_=1., log_normal=True)