gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
statistics.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
19 
20 #include <cmath>
21 #include <limits>
22 #include <fstream>
23 #include <iomanip>
24 
29 
30 namespace Gambit
31 {
32 
33  namespace Stats
34  {
35 
37  const double logmin = log(std::numeric_limits<double>::min());
38 
41  {
42  const int npts = 1001;
43  const int n_theory_errs = 5;
44  const double obs = 10;
45  const double obserr = 1.0;
46  const double base_theory_err = 0.25;
47  const double obsrelerr = obserr/obs;
48  const double pred_min = 0;
49  const double pred_max = 20;
50  double theory_err[n_theory_errs] = {0, 1, 2, 4, 8};
51  double theory_relerr[n_theory_errs];
52  for (int j = 0; j < n_theory_errs; ++j)
53  {
54  theory_err[j] *= base_theory_err;
55  theory_relerr[j] = theory_err[j]/obs;
56  }
57 
58  std::ofstream myfile;
59  myfile.open ("stats_tests.txt");
60  myfile << "# Theory errors:";
61  for (int j = 0; j < n_theory_errs; ++j) myfile << " " << theory_err[j];
62  myfile << endl << "# Fractional (relative) theory errors:";
63  for (int j = 0; j < n_theory_errs; ++j) myfile << " " << theory_relerr[j];
64  myfile << endl << "# mu ";
65  for (int j = 0; j < n_theory_errs; ++j)
66  {
67  myfile << "ln_Gprof[" << j << "] "
68  << "ln_Gmarg[" << j << "] "
69  << "ln_Gupprof[" << j << "] "
70  << "ln_Gupmarg[" << j << "] "
71  << "ln_Gdownprof[" << j << "] "
72  << "ln_Gdownmarg[" << j << "] "
73  << "ln_LNprof[" << j << "] "
74  << "ln_LNmarg[" << j << "] "
75  << "ln_LNprof_rel[" << j << "]"
76  << "ln_LNmarg_rel[" << j << "]"
77  << "Gprof[" << j << "] "
78  << "Gmarg[" << j << "] "
79  << "Gupprof[" << j << "] "
80  << "Gupmarg[" << j << "] "
81  << "Gdownprof[" << j << "] "
82  << "Gdownmarg[" << j << "] "
83  << "LNprof[" << j << "] "
84  << "LNmarg[" << j << "] "
85  << "LNprof_rel[" << j << "] "
86  << "LNmarg_rel[" << j << "] ";
87  }
88  myfile << endl << std::setprecision(8) << std::scientific;
89  for (int i = 0; i < npts; ++i)
90  {
91  double pred = pred_min + double(i)/double(npts-1)*(pred_max - pred_min);
92  myfile << pred << " ";
93  for (int j = 0; j < n_theory_errs; ++j)
94  {
95  double Gprof = gaussian_loglikelihood(pred, obs, theory_err[j], obserr, true);
96  double Gmarg = gaussian_loglikelihood(pred, obs, theory_err[j], obserr, false);
97  double Gupprof = gaussian_upper_limit(pred, obs, theory_err[j], obserr, true);
98  double Gupmarg = gaussian_upper_limit(pred, obs, theory_err[j], obserr, false);
99  double Gdownprof = gaussian_lower_limit(pred, obs, theory_err[j], obserr, true);
100  double Gdownmarg = gaussian_lower_limit(pred, obs, theory_err[j], obserr, false);
101  double LNprof = lognormal_loglikelihood(pred, obs, theory_err[j], obserr, true);
102  double LNmarg = lognormal_loglikelihood(pred, obs, theory_err[j], obserr, false);
103  double LNprof_rel = lognormal_loglikelihood_relerror(pred, obs, theory_relerr[j], obsrelerr, true);
104  double LNmarg_rel = lognormal_loglikelihood_relerror(pred, obs, theory_relerr[j], obsrelerr, false);
105  myfile << Gprof << " "
106  << Gmarg << " "
107  << Gupprof << " "
108  << Gupmarg << " "
109  << Gdownprof << " "
110  << Gdownmarg << " "
111  << LNprof << " "
112  << LNmarg << " "
113  << LNprof_rel << " "
114  << LNmarg_rel << " "
115  << exp(Gprof) << " "
116  << exp(Gmarg) << " "
117  << exp(Gupprof) << " "
118  << exp(Gupmarg) << " "
119  << exp(Gdownprof) << " "
120  << exp(Gdownmarg) << " "
121  << exp(LNprof) << " "
122  << exp(LNmarg) << " "
123  << exp(LNprof_rel) << " "
124  << exp(LNmarg_rel) << " ";
125  }
126  myfile << endl;
127  }
128  myfile.close();
129  exit(1);
130  }
131 
133  double gaussian_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
134  {
135  //Debug code (used to produce plots as per Section 8 of the main GAMBIT paper).
136  //static bool first = true;
137  //if (first) {first = false; test_likelihoods();}
138 
139  if (theoryerr < 0) utils_error().raise(LOCAL_INFO, "Theory uncertainty cannot be negative.");
140  if (obserr <= 0) utils_error().raise(LOCAL_INFO, "Observational uncertainty must be non-zero and positive.");
141 
142  // If no theory error is given, go back to the standard chi^2 without systematics.
143  if (theoryerr == 0) profile_systematics = false;
144 
145  double errsq = theoryerr*theoryerr + obserr*obserr;
146  double chi2 = -0.5*pow(theory-obs,2)/errsq;
147  double norm = profile_systematics ? log(2.0*pi*obserr*theoryerr) : 0.5*log(2.0*pi*errsq);
148 
149  return chi2 - norm;
150  }
151 
154  double lognormal_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
155  {
156  return lognormal_loglikelihood_relerror(theory, obs, theoryerr/theory, obserr/obs, profile_systematics);
157  }
158 
161  double lognormal_loglikelihood_relerror(double theory, double obs, double reltheoryerr, double relobserr, bool profile_systematics)
162  {
163  if (reltheoryerr < 0) utils_error().raise(LOCAL_INFO, "Theory uncertainty cannot be negative.");
164  if (relobserr <= 0) utils_error().raise(LOCAL_INFO, "Observational uncertainty must be non-zero and positive.");
165 
166  // Don't allow observed or predicted values of less than or equal to 0
167  if (obs <= 1e2*std::numeric_limits<double>::min()) return std::numeric_limits<double>::lowest();
168  if (theory <= 1e2*std::numeric_limits<double>::min()) return std::numeric_limits<double>::lowest();
169 
170  // If no theory error is given, go back to the standard log-normal likelihood without systematics.
171  if (reltheoryerr == 0) profile_systematics = false;
172 
173  double log_obs_on_theory = log(obs / theory);
174  double obserr_prime = log1p(relobserr);
175  double theoryerr_prime = log1p(reltheoryerr);
176  double errsq = obserr_prime*obserr_prime + theoryerr_prime*theoryerr_prime;
177  double chi2 = -log(obs) - 0.5*pow(log_obs_on_theory,2)/errsq;
178  double norm = profile_systematics ? log(2.0*pi*obserr_prime*theoryerr_prime) : 0.5*log(2.0*pi*errsq);
179  double extra = profile_systematics ? theoryerr_prime*theoryerr_prime*(0.5*obserr_prime*obserr_prime - log_obs_on_theory)/errsq : 0.0;
180  return chi2 + extra - norm;
181  }
182 
184  double gaussian_upper_limit(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
185  {
186  if (theoryerr < 0) utils_error().raise(LOCAL_INFO, "Theory uncertainty cannot be negative.");
187  if (obserr <= 0) utils_error().raise(LOCAL_INFO, "Observational uncertainty must be non-zero and positive.");
188 
189  double loglike;
190  double errsq = theoryerr*theoryerr + obserr*obserr;
191 
192  // Deal with the special case that no theory error has actually been given. Revert to basic limit likelihood.
193  if (theoryerr == 0)
194  {
195  double prefactor = -0.5*log(2.0*pi*errsq);
196  return prefactor - (theory > obs ? 0.5*pow(theory-obs,2)/errsq : 0.0);
197  }
198 
199  // Switch according to whether or not we are profiling.
200  if (profile_systematics)
201  {
202  double prefactor = -log(2.0*pi*theoryerr*obserr);
203  loglike = prefactor - (theory > obs? 0.5*pow(theory-obs,2)/errsq : 0.0);
204  }
205  else
206  {
207  double prefactor = -0.5*log(8.0*pi);
208  double diff = obs - theory;
209  double like = exp(-0.5*pow(diff,2)/errsq)/sqrt(errsq);
210  like *= erfc(obserr * diff / (theoryerr * sqrt(2.0*errsq)));
211  like += erfc(-diff/(sqrt(2.0)*theoryerr))/obserr;
212  if (like < 0) utils_error().raise(LOCAL_INFO, "Marginalised Gaussian limit likelihood went negative!");
213  loglike = prefactor + (like == 0 ? logmin : log(like));
214  }
215 
216  return loglike;
217 
218  }
219 
221  double gaussian_lower_limit(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
222  {
223  return gaussian_upper_limit(-theory, -obs, theoryerr, obserr, profile_systematics);
224  }
225 
226  }
227 
228 }
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
EXPORT_SYMBOLS error & utils_error()
Utility errors.
LOCAL_INFO macro.
#define LOCAL_INFO
Definition: local_info.hpp:34
double gaussian_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is Gaussian...
Definition: statistics.cpp:133
void test_likelihoods()
A simple tester routine for the likelihood functions.
Definition: statistics.cpp:40
double lognormal_loglikelihood_relerror(double theory, double obs, double reltheoryerr, double relobserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is log-norm...
Definition: statistics.cpp:161
Declarations of statistical utilities.
const double logmin
Minimum finite result returnable from log(double x);.
Definition: statistics.cpp:37
const double pi
double gaussian_upper_limit(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a gaussian log-likelihood for an upper limit.
Definition: statistics.cpp:184
double gaussian_lower_limit(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a gaussian log-likelihood for a lower limit.
Definition: statistics.cpp:221
Basic set of known mathematical and physical constants for GAMBIT.
double pow(const double &a)
Outputs a^i.
Exception objects required for standalone compilation.
double lognormal_loglikelihood(double theory, double obs, double theoryerr, double obserr, bool profile_systematics)
Use a detection to compute a simple chi-square-like log likelihood, for the case when obs is log-norm...
Definition: statistics.cpp:154
TODO: see if we can use this one:
Definition: Analysis.hpp:33