conditionalgaussian.cpp
Go to the documentation of this file.
00001 // $Id: conditionalgaussian.cpp 29890 2009-02-02 10:22:01Z tdelaet $
00002 // Copyright (C) 2003 Klaas Gadeyne <first dot last at gmail dot com>
00003 // Copyright (C) 2008 Tinne De Laet <first dot last at mech dot kuleuven dot be>
00004 //
00005 // This program is free software; you can redistribute it and/or modify
00006 // it under the terms of the GNU Lesser General Public License as published by
00007 // the Free Software Foundation; either version 2.1 of the License, or
00008 // (at your option) any later version.
00009 //
00010 // This program is distributed in the hope that it will be useful,
00011 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013 // GNU Lesser General Public License for more details.
00014 //
00015 // You should have received a copy of the GNU Lesser General Public License
00016 // along with this program; if not, write to the Free Software
00017 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00018 //
00019 
00020 #include "conditionalgaussian.h"
00021 #include <cmath>
00022 #include "../wrappers/rng/rng.h"
00023 
00024 namespace BFL
00025 {
00026   using namespace MatrixWrapper;
00027 
00028   ConditionalGaussian::ConditionalGaussian(int dim,
00029                                            int num_conditional_arguments)
00030     : ConditionalPdf<ColumnVector,ColumnVector>(dim, num_conditional_arguments)
00031     , _diff(dim)
00032     , _Mu(dim)
00033     , _Low_triangle(dim,dim)
00034     , _samples(dim)
00035     , _SampleValue(dim)
00036   {}
00037 
00039   ConditionalGaussian::~ConditionalGaussian(){}
00040 
00041   //Clone function
00042   ConditionalGaussian* ConditionalGaussian::Clone() const
00043   {
00044       return new ConditionalGaussian(*this);
00045   }
00046 
00047   Probability
00048   ConditionalGaussian::ProbabilityGet(const ColumnVector& input) const
00049   {
00050     // Update Mu
00051     _Mu = ExpectedValueGet();
00052     _diff = input - _Mu;
00053 
00054     Probability temp = _diff.transpose() * (ColumnVector)(CovarianceGet().inverse() * _diff);
00055     Probability result = exp(-0.5 * temp) / sqrt(pow(M_PI*2,(double)DimensionGet()) * CovarianceGet().determinant());
00056     return result;
00057   }
00058 
00059   bool
00060   ConditionalGaussian::SampleFrom (vector<Sample<ColumnVector> >& samples, const int num_samples, int method, void * args) const
00061   {
00062     return Pdf<ColumnVector>::SampleFrom(samples, num_samples, method, args);
00063   }
00064 
00065   bool
00066   ConditionalGaussian::SampleFrom (Sample<ColumnVector>& sample, int method, void * args) const
00067   {
00068     // Sampling from a Gaussian is simple if DIMENSION = 1 or 2 (and the
00069     // 2 variables are independant!)
00070     // Then we can use inversion sampling (Box-Muller method)
00071     // So for 1D, we use Box-Muller, else we use the cholesky method
00072     // These are both methods that don't require any arguments
00073 
00074     // Update mu
00075     _Mu = ExpectedValueGet();
00076 
00077     switch(method)
00078       {
00079       case DEFAULT: // Cholesky, see althere (bad implementation)
00080         {
00081           bool result = CovarianceGet().cholesky_semidefinite(_Low_triangle);
00082           for (unsigned int j=1; j < DimensionGet()+1; j++){_samples(j) = rnorm(0,1);}
00083           _SampleValue = (_Low_triangle * _samples) + _Mu;
00084           sample.ValueSet(_SampleValue);
00085           return result;
00086         }
00087       case BOXMULLER: 
00088         {
00089           cerr << "Box-Muller not implemented yet!" << endl;
00090           return false;
00091         }
00092       case CHOLESKY: // Cholesky Sampling
00093         {
00094           bool result = CovarianceGet().cholesky_semidefinite(_Low_triangle);
00095           /* For now we keep it simple, and use the scythe library
00096              (although wrapped) with the uRNG that it uses itself only */
00097           /* Sample Gaussian._dimension samples from univariate
00098              gaussian This could be done using several available
00099              libraries, combined with different uniform RNG.  Both the
00100              library to be used and the uRNG could be implemented as
00101              #ifdef conditions, although I'm sure there must exist a
00102              cleaner way to implement this!
00103           */
00104           for (unsigned int j=1; j < DimensionGet()+1; j++) _samples(j) = rnorm(0,1);
00105           _SampleValue = (_Low_triangle * _samples) + _Mu;
00106           sample.ValueSet(_SampleValue);
00107           return result;
00108         }
00109       default:
00110         cerr << "Conditional Gaussian: Sampling method " << method
00111              << "not implemented yet!" << endl;
00112         return false;
00113       }
00114   }
00115 
00116 } // End namespace


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 Sun Oct 5 2014 22:29:52