gaussian.cpp
Go to the documentation of this file.
1 // $Id$
2 // Copyright (C) 2002 Klaas Gadeyne <first dot last at gmail dot com>
3 // Copyright (C) 2008 Tinne De Laet <first dot last at mech dot kuleuven dot be>
4 //
5 // This program is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published by
7 // the Free Software Foundation; either version 2.1 of the License, or
8 // (at your option) any later version.
9 //
10 // This program is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with this program; if not, write to the Free Software
17 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18 //
19 #include "gaussian.h"
20 
21 #include "../wrappers/rng/rng.h" // Wrapper around several rng libraries
22 
23 #include <cmath> // For sqrt and exp
24 #include <cassert>
25 
26 namespace BFL
27 {
28  using namespace MatrixWrapper;
29 
30  Gaussian::Gaussian (const ColumnVector& m, const SymmetricMatrix& s)
31  : Pdf<ColumnVector> ( m.rows() )
32  , _diff(DimensionGet())
33  , _tempColumn(DimensionGet())
34  , _samples(DimensionGet())
35  , _sampleValue(DimensionGet())
36  , _Low_triangle(DimensionGet(),DimensionGet())
37  {
38  // check if inputs are consistent
39  assert (m.rows() == s.columns());
40  _Mu = m;
41  _Sigma = s;
42  _Sigma_inverse.resize(DimensionGet());
43  _Sigma_changed = true;
44  }
45 
46  Gaussian::Gaussian (int dimension)
47  : Pdf<ColumnVector>(dimension)
48  , _diff(dimension)
49  , _tempColumn(DimensionGet())
50  , _samples(dimension)
51  , _sampleValue(dimension)
52  , _Low_triangle(dimension,dimension)
53  {
54  _Mu.resize(dimension);
55  _Sigma.resize(dimension);
56  _Sigma_inverse.resize(dimension);
57  _Sigma_changed = true;
58  }
59 
61 
62  std::ostream& operator<< (std::ostream& os, const Gaussian& g)
63  {
64  os << "\nMu:\n" << g.ExpectedValueGet()
65  << "\nSigma:\n" << g.CovarianceGet() << endl;
66  return os;
67  }
68 
69  //Clone function
71  {
72  return new Gaussian(*this);
73  }
74 
75  Probability Gaussian::ProbabilityGet(const ColumnVector& input) const
76  {
77  // only calculate these variables if sigma has changed
78  if (_Sigma_changed){
79  _Sigma_changed = false;
80  _Sigma_inverse = _Sigma.inverse();
81  _sqrt_pow = 1 / sqrt(pow(M_PI*2,(double)DimensionGet()) * _Sigma.determinant());
82  }
83 
84  _diff = input;
85  _diff -= _Mu;
86  _Sigma_inverse.multiply(_diff,_tempColumn);
87  //_tempColumn = _Sigma_inverse * _diff;
88  Probability temp = _diff.transpose() * _tempColumn;
89  //Probability temp = _diff.transpose() * (_Sigma_inverse * _diff);
90  Probability result = exp(-0.5 * temp) * _sqrt_pow;
91  return result;
92  }
93 
94  // Redefined for optimal performance. Eg. do Cholesky decomposition
95  // only once when drawing multiple samples at once!
96  // See method below for more info regarding the algorithms
97  bool
98  Gaussian::SampleFrom (vector<Sample<ColumnVector> >& list_samples, const int num_samples, int method, void * args) const
99  {
100  list_samples.resize(num_samples); // will break real-timeness if list_samples.size()!=num_samples
101  vector<Sample<ColumnVector> >::iterator rit = list_samples.begin();
102  switch(method)
103  {
104  case DEFAULT: // Cholesky Sampling
105  case CHOLESKY:
106  {
107  bool result = _Sigma.cholesky_semidefinite(_Low_triangle);
108  while (rit != list_samples.end())
109  {
110  for (unsigned int j=1; j < DimensionGet()+1; j++) _samples(j) = rnorm(0,1);
112  _sampleValue += this->_Mu;
113  rit->ValueSet(_sampleValue);
114  rit++;
115  }
116  return result;
117  }
118  case BOXMULLER: // Implement box-muller here
119  // Only for univariate distributions.
120  return false;
121  default:
122  return false;
123  }
124  }
125 
126 
127  bool
128  Gaussian::SampleFrom (Sample<ColumnVector>& one_sample, int method, void * args) const
129  {
130  /* Exact i.i.d. samples from a Gaussian can be drawn in several
131  ways:
132  - if the DimensionGet() = 1 or 2 (and the 2 variables are
133  independant), we can use inversion sampling (Box-Muller
134  method)
135  - For larger dimensions, we use can use the Cholesky method or
136  an approached based on conditional distributions.
137  (see ripley87, P.98 (bibtex below)). The Cholesky method is
138  generally preferred and the only one implemented for now.
139  */
140  switch(method)
141  {
142  case DEFAULT: // Cholesky Sampling, see eg.
143  case CHOLESKY: // Cholesky Sampling, see eg.
144  /*
145  @Book{ ripley87,
146  author = {Ripley, Brian D.},
147  title = {Stochastic Simulation},
148  publisher = {John Wiley and Sons},
149  year = {1987},
150  annote = {ISBN 0271-6356, WBIB 1 519.245}
151  }
152  p.98
153  */
154  {
155  bool result = _Sigma.cholesky_semidefinite(_Low_triangle);
156 
157  /* For now we keep it simple, and use the scythe library
158  (although wrapped) with the uRNG that it uses itself only */
159  /* Sample Gaussian._dimension samples from univariate
160  gaussian This could be done using several available
161  libraries, combined with different uniform RNG. Both the
162  library to be used and the uRNG could be implemented as
163  #ifdef conditions, although I'm sure there must exist a
164  cleaner way to implement this!
165  */
166  for (unsigned int j=1; j < DimensionGet()+1; j++) _samples(j) = rnorm(0,1);
167  _sampleValue = (_Low_triangle * _samples) + this->_Mu;
168  one_sample.ValueSet(_sampleValue);
169  return result;
170  }
171  case BOXMULLER: // Implement box-muller here
172  // Only for univariate distributions.
173  return false;
174  default:
175  return false;
176  }
177  }
178 
179 
180  ColumnVector
182  {
183  return _Mu;
184  }
185 
186  SymmetricMatrix
188  {
189  return _Sigma;
190  }
191 
192  void
193  Gaussian::ExpectedValueSet (const ColumnVector& mu)
194  {
195  _Mu = mu;
196  if (this->DimensionGet() == 0)
197  {
198  this->DimensionSet(mu.rows());
199  }
200  assert(this->DimensionGet() == mu.rows());
201  }
202 
203  void
204  Gaussian::CovarianceSet (const SymmetricMatrix& cov)
205  {
206  _Sigma = cov;
207  _Sigma_changed = true;
208  if (this->DimensionGet() == 0)
209  {
210  this->DimensionSet(cov.rows());
211  }
212  assert(this->DimensionGet() == cov.rows());
213  }
214 
215 } // End namespace BFL
void ExpectedValueSet(const MatrixWrapper::ColumnVector &mu)
Set the Expected Value.
Definition: gaussian.cpp:193
Class PDF: Virtual Base class representing Probability Density Functions.
Definition: pdf.h:53
virtual Probability ProbabilityGet(const MatrixWrapper::ColumnVector &input) const
Get the probability of a certain argument.
Definition: gaussian.cpp:75
#define BOXMULLER
ColumnVector _samples
Definition: gaussian.h:40
#define DEFAULT
friend std::ostream & operator<<(std::ostream &os, const Gaussian &g)
output stream for Gaussian
Definition: gaussian.cpp:62
MatrixWrapper::ColumnVector _Mu
Definition: gaussian.h:30
double rnorm(const double &mu, const double &sigma)
Class representing Gaussian (or normal density)
Definition: gaussian.h:27
bool _Sigma_changed
Definition: gaussian.h:34
MatrixWrapper::SymmetricMatrix _Sigma
Definition: gaussian.h:31
void ValueSet(const T &value)
Set the value of the Sample.
MatrixWrapper::SymmetricMatrix _Sigma_inverse
Definition: gaussian.h:35
virtual MatrixWrapper::SymmetricMatrix CovarianceGet() const
Get the Covariance Matrix E[(x - E[x])^2] of the Analytic pdf.
Definition: gaussian.cpp:187
bool SampleFrom(vector< Sample< MatrixWrapper::ColumnVector > > &list_samples, const int num_samples, int method=DEFAULT, void *args=NULL) const
virtual ~Gaussian()
Default Copy Constructor will do.
Definition: gaussian.cpp:60
Gaussian(const MatrixWrapper::ColumnVector &Mu, const MatrixWrapper::SymmetricMatrix &Sigma)
Constructor.
#define M_PI
Definition: bfl_constants.h:26
Matrix _Low_triangle
Definition: gaussian.h:42
unsigned int DimensionGet() const
Get the dimension of the argument.
virtual Gaussian * Clone() const
Clone function.
Definition: gaussian.cpp:70
ColumnVector _diff
Definition: gaussian.h:37
#define CHOLESKY
ColumnVector _tempColumn
Definition: gaussian.h:38
void CovarianceSet(const MatrixWrapper::SymmetricMatrix &cov)
Set the Covariance Matrix.
Definition: gaussian.cpp:204
ColumnVector _sampleValue
Definition: gaussian.h:41
double _sqrt_pow
Definition: gaussian.h:36
Class representing a probability (a double between 0 and 1)
Definition: bfl_constants.h:39
virtual void DimensionSet(unsigned int dim)
Set the dimension of the argument.
virtual MatrixWrapper::ColumnVector ExpectedValueGet() const
Get the expected value E[x] of the pdf.
Definition: gaussian.cpp:181


bfl
Author(s): Klaas Gadeyne, Wim Meeussen, Tinne Delaet and many others. See web page for a full contributor list. ROS package maintained by Wim Meeussen.
autogenerated on Mon Feb 28 2022 21:56:33