gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
gaussian.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
24 
25 #ifndef __PRIOR_GAUSSIAN_HPP__
26 #define __PRIOR_GAUSSIAN_HPP__
27 
28 #include <algorithm>
29 #include <cmath>
30 #include <string>
31 #include <unordered_map>
32 #include <vector>
33 
37 
38 #include <boost/math/special_functions/erf.hpp>
39 
40 namespace Gambit
41 {
42  namespace Priors
43  {
52  class Gaussian : public BasePrior
53  {
54  private:
55  std::vector <double> mu;
56  mutable Cholesky col;
57 
58  public:
59  // Constructor defined in gaussian.cpp
60  Gaussian(const std::vector<std::string>&, const Options&);
61 
63  void transform(const std::vector <double> &unitpars, std::unordered_map<std::string, double> &outputMap) const
64  {
65  std::vector<double> vec(unitpars.size());
66 
67  auto v_it = vec.begin();
68  for (auto elem_it = unitpars.begin(), elem_end = unitpars.end(); elem_it != elem_end; elem_it++, v_it++)
69  {
70  *v_it = M_SQRT2 * boost::math::erf_inv(2. * (*elem_it) - 1.);
71  }
72 
73  col.ElMult(vec);
74 
75  v_it = vec.begin();
76  auto m_it = mu.begin();
77  for (auto str_it = param_names.begin(), str_end = param_names.end(); str_it != str_end; str_it++)
78  {
79  outputMap[*str_it] = *(v_it++) + *(m_it++);
80  }
81  }
82 
83  std::vector<double> inverse_transform(const std::unordered_map<std::string, double> &physical) const override
84  {
85  // subtract mean
86  std::vector<double> central;
87  for (int i = 0, n = this->size(); i < n; i++)
88  {
89  central.push_back(physical.at(param_names[i]) - mu[i]);
90  }
91 
92  // invert rotation by Cholesky matrix
93  std::vector<double> rotated = col.invElMult(central);
94 
95  // now diagonal; invert Gaussian CDF
96  std::vector<double> u;
97  for (const auto& v : rotated)
98  {
99  u.push_back(0.5 * (boost::math::erf(v / M_SQRT2) + 1.));
100  }
101  return u;
102  }
103 
104  double operator()(const std::vector<double> &vec) const
105  {
106  static double norm = 0.5 * std::log(2. * M_PI * std::pow(col.DetSqrt(), 2));
107  return -0.5 * col.Square(vec, mu) - norm;
108  }
109  };
110 
111  LOAD_PRIOR(gaussian, Gaussian)
112 
113  } // namespace Priors
114 } // namespace Gambit
115 
116 #endif // __PRIOR_GAUSSIAN_HPP__
Abstract base class for priors.
Definition: base_prior.hpp:40
declaration for scanner module
std::vector< double > mu
Definition: gaussian.hpp:55
double operator()(const std::vector< double > &vec) const
Log of PDF density.
Definition: gaussian.hpp:104
std::vector< std::string > param_names
Definition: base_prior.hpp:46
double Square(const std::vector< double > &y, const std::vector< double > &y0)
Definition: cholesky.hpp:106
Prior object construction routines.
Declarations for the YAML options class.
unsigned int size() const
Definition: base_prior.hpp:77
Multi-dimensional Gaussian prior.
Definition: gaussian.hpp:52
double DetSqrt()
Definition: cholesky.hpp:127
std::vector< double > inverse_transform(const std::unordered_map< std::string, double > &physical) const override
Transform from parameter back to unit hypercube.
Definition: gaussian.hpp:83
std::vector< double > invElMult(const std::vector< double > &y) const
x = L^-1 y where L is the lower-diagonal Cholesky matrix
Definition: cholesky.hpp:91
void transform(const std::vector< double > &unitpars, std::unordered_map< std::string, double > &outputMap) const
Transformation from unit interval to the Gaussian.
Definition: gaussian.hpp:63
RangePrior1D< flatprior > LOAD_PRIOR(cos, RangePrior1D< cosprior >) LOAD_PRIOR(sin
double pow(const double &a)
Outputs a^i.
void ElMult(std::vector< double > &y) const
Definition: cholesky.hpp:69
std::vector< T > vec(std::vector< T > vector)
Definition: daFunk.hpp:142
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.
Gaussian(const std::vector< std::string > &, const Options &)
Definition: gaussian.cpp:31