conditionalgaussian.cpp
Go to the documentation of this file.
1 // \$Id\$
2 // Copyright (C) 2003 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
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
20 #include "conditionalgaussian.h"
21 #include <cmath>
22 #include "../wrappers/rng/rng.h"
23
24 namespace BFL
25 {
26  using namespace MatrixWrapper;
27
29  int num_conditional_arguments)
30  : ConditionalPdf<ColumnVector,ColumnVector>(dim, num_conditional_arguments)
31  , _diff(dim)
32  , _Mu(dim)
33  , _Low_triangle(dim,dim)
34  , _samples(dim)
35  , _SampleValue(dim)
36  {}
37
40
41  //Clone function
43  {
44  return new ConditionalGaussian(*this);
45  }
46
48  ConditionalGaussian::ProbabilityGet(const ColumnVector& input) const
49  {
50  // Update Mu
52  _diff = input - _Mu;
53
54  Probability temp = _diff.transpose() * (ColumnVector)(CovarianceGet().inverse() * _diff);
55  Probability result = exp(-0.5 * temp) / sqrt(pow(M_PI*2,(double)DimensionGet()) * CovarianceGet().determinant());
56  return result;
57  }
58
59  bool
60  ConditionalGaussian::SampleFrom (vector<Sample<ColumnVector> >& samples, const int num_samples, int method, void * args) const
61  {
62  return Pdf<ColumnVector>::SampleFrom(samples, num_samples, method, args);
63  }
64
65  bool
66  ConditionalGaussian::SampleFrom (Sample<ColumnVector>& sample, int method, void * args) const
67  {
68  // Sampling from a Gaussian is simple if DIMENSION = 1 or 2 (and the
69  // 2 variables are independant!)
70  // Then we can use inversion sampling (Box-Muller method)
71  // So for 1D, we use Box-Muller, else we use the cholesky method
72  // These are both methods that don't require any arguments
73
74  // Update mu
76
77  switch(method)
78  {
79  case DEFAULT: // Cholesky, see althere (bad implementation)
80  {
81  bool result = CovarianceGet().cholesky_semidefinite(_Low_triangle);
82  for (unsigned int j=1; j < DimensionGet()+1; j++){_samples(j) = rnorm(0,1);}
84  sample.ValueSet(_SampleValue);
85  return result;
86  }
87  case BOXMULLER:
88  {
89  cerr << "Box-Muller not implemented yet!" << endl;
90  return false;
91  }
92  case CHOLESKY: // Cholesky Sampling
93  {
94  bool result = CovarianceGet().cholesky_semidefinite(_Low_triangle);
95  /* For now we keep it simple, and use the scythe library
96  (although wrapped) with the uRNG that it uses itself only */
97  /* Sample Gaussian._dimension samples from univariate
98  gaussian This could be done using several available
99  libraries, combined with different uniform RNG. Both the
100  library to be used and the uRNG could be implemented as
101  #ifdef conditions, although I'm sure there must exist a
102  cleaner way to implement this!
103  */
104  for (unsigned int j=1; j < DimensionGet()+1; j++) _samples(j) = rnorm(0,1);
106  sample.ValueSet(_SampleValue);
107  return result;
108  }
109  default:
110  cerr << "Conditional Gaussian: Sampling method " << method
111  << "not implemented yet!" << endl;
112  return false;
113  }
114  }
115
116 } // End namespace
#define BOXMULLER
#define DEFAULT
double rnorm(const double &mu, const double &sigma)
virtual Probability ProbabilityGet(const MatrixWrapper::ColumnVector &input) const
Get the probability of a certain argument.
void ValueSet(const T &value)
Set the value of the Sample.
virtual bool SampleFrom(Sample< MatrixWrapper::ColumnVector > &sample, int method=DEFAULT, void *args=NULL) const
virtual ConditionalGaussian * Clone() const
Clone function.
Abstract Class representing conditional Pdfs P(x | ...)
virtual bool SampleFrom(vector< Sample< T > > &list_samples, const unsigned int num_samples, int method=DEFAULT, void *args=NULL) const
Draw multiple samples from the Pdf (overloaded)
#define M_PI
Definition: bfl_constants.h:26
#define CHOLESKY
virtual ~ConditionalGaussian()
Destructor.
ConditionalGaussian(int dim=0, int num_conditional_arguments=0)
Constructor.
Class representing a probability (a double between 0 and 1)
Definition: bfl_constants.h:39
unsigned int DimensionGet() const
Get the dimension of the argument.
Abstract Class representing all Conditional gaussians.
virtual MatrixWrapper::SymmetricMatrix CovarianceGet() const
Get the Covariance Matrix E[(x - E[x])^2] of the Analytic pdf.
virtual MatrixWrapper::ColumnVector ExpectedValueGet() const
Get the expected value E[x] of the pdf.

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 Jun 10 2019 12:47:58