EigenMultiVariateNormal.hpp
Go to the documentation of this file.
00001 #ifndef __EIGENMULTIVARIATENORMAL_HPP
00002 #define __EIGENMULTIVARIATENORMAL_HPP
00003 
00004 //found at http://lost-found-wandering.blogspot.de/2011/05/sampling-from-multivariate-normal-in-c.html
00005 
00006 #include <Eigen/Dense>
00007 
00008 #include <math.h>
00009 
00010 #include <boost/random/mersenne_twister.hpp>
00011 #include <boost/random/normal_distribution.hpp>
00012 #include <boost/random/variate_generator.hpp>
00013 
00020 template<typename _Scalar, int _size>
00021 class EigenMultivariateNormal
00022 {
00023     boost::mt19937 rng;    // The uniform pseudo-random algorithm
00024     boost::normal_distribution<_Scalar> norm;  // The gaussian combinator
00025     boost::variate_generator<boost::mt19937&,boost::normal_distribution<_Scalar> >
00026        randN; // The 0-mean unit-variance normal generator
00027 
00028     Eigen::Matrix<_Scalar,_size,_size> rot;
00029     Eigen::Matrix<_Scalar,_size,1> scl;
00030 
00031     Eigen::Matrix<_Scalar,_size,1> mean;
00032 
00033 public:
00034     EigenMultivariateNormal(const Eigen::Matrix<_Scalar,_size,1>& meanVec,
00035         const Eigen::Matrix<_Scalar,_size,_size>& covarMat)
00036         : randN(rng,norm)
00037     {
00038         setCovar(covarMat);
00039         setMean(meanVec);
00040     }
00041 
00042     void setCovar(const Eigen::Matrix<_Scalar,_size,_size>& covarMat)
00043     {
00044         Eigen::SelfAdjointEigenSolver<Eigen::Matrix<_Scalar,_size,_size> >
00045            eigenSolver(covarMat);
00046         rot = eigenSolver.eigenvectors();
00047         scl = eigenSolver.eigenvalues();
00048         for (int ii=0;ii<_size;++ii) {
00049             scl(ii,0) = sqrt(scl(ii,0));
00050         }
00051     }
00052 
00053     void setMean(const Eigen::Matrix<_Scalar,_size,1>& meanVec)
00054     {
00055         mean = meanVec;
00056     }
00057 
00058     void nextSample(Eigen::Matrix<_Scalar,_size,1>& sampleVec)
00059     {
00060         for (int ii=0;ii<_size;++ii) {
00061             sampleVec(ii,0) = randN()*scl(ii,0);
00062         }
00063         sampleVec = rot*sampleVec + mean;
00064     }
00065 
00066 };
00067 
00068 
00069 
00073 template<typename _Scalar>
00074 class UnivariateNormal
00075 {
00076     boost::mt19937 rng;    // The uniform pseudo-random algorithm
00077     boost::normal_distribution<_Scalar> norm;  // The gaussian combinator
00078     boost::variate_generator<boost::mt19937&,boost::normal_distribution<_Scalar> >
00079        randN; // The 0-mean unit-variance normal generator
00080 
00081     _Scalar mean_;
00082     _Scalar stdev_;
00083 
00084 public:
00085 
00086     UnivariateNormal(const _Scalar & mean, const _Scalar &var)
00087         : randN(rng,norm)
00088     {
00089         setVar(var);
00090         setMean(mean);
00091     }
00092 
00093     void setVar(const _Scalar & var)
00094     {
00095         stdev_ = sqrt(fabs(var));
00096     }
00097 
00098     void setMean(const _Scalar & mean)
00099     {
00100         mean_ = mean;
00101     }
00102 
00103     void nextSample(_Scalar &sample)
00104     {
00105         sample = randN() * stdev_ + mean_;
00106     }
00107 
00108 };
00109 
00110 
00111 
00112 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Properties Friends


uncertain_tf
Author(s): Thomas Ruehr
autogenerated on Thu May 23 2013 17:49:36