[3] | 1 | #ifndef __DISTRIBUTIONS_HPP__
|
---|
| 2 | #define __DISTRIBUTIONS_HPP__
|
---|
| 3 |
|
---|
| 4 | #include "ublas.hpp"
|
---|
| 5 | #include <boost/math/distributions.hpp>
|
---|
| 6 |
|
---|
| 7 | namespace math{
|
---|
| 8 | namespace distributions{
|
---|
| 9 |
|
---|
| 10 |
|
---|
| 11 | template <class RealType>
|
---|
| 12 | inline RealType pdf(const boost::math::normal_distribution<RealType> & dist , boost::numeric::ublas::vector<RealType> & v){
|
---|
| 13 | return pdf(dist,v[0]);
|
---|
| 14 | }
|
---|
| 15 |
|
---|
| 16 | template <class RealType>
|
---|
| 17 | inline RealType pdf(const boost::math::uniform_distribution<RealType> & dist , boost::numeric::ublas::vector<RealType> & v){
|
---|
| 18 | return pdf(dist,v[0]);
|
---|
| 19 | }
|
---|
| 20 |
|
---|
| 21 |
|
---|
| 22 | /*!
|
---|
| 23 | *\class multivariate_normal_distribution
|
---|
| 24 | *\brief This clas describes a multivariate normal distribution
|
---|
| 25 | */
|
---|
| 26 | template < class RealType =double >
|
---|
| 27 | class multivariate_normal_distribution{
|
---|
| 28 | public :
|
---|
| 29 |
|
---|
| 30 | /*!
|
---|
| 31 | *\brief Default constructor
|
---|
| 32 | */
|
---|
| 33 | multivariate_normal_distribution(boost::numeric::ublas::vector<RealType> mean,boost::numeric::ublas::matrix<RealType> cov){
|
---|
| 34 |
|
---|
| 35 | using namespace boost::numeric::ublas;
|
---|
| 36 |
|
---|
| 37 | if(mean.size()==cov.size1() && mean.size()==cov.size2()){
|
---|
| 38 | m_mean=mean;
|
---|
| 39 | m_cov=cov;
|
---|
| 40 | }else throw math_error("Multivariate normal distribution : the mean vector of covariance matrix must have the same size");
|
---|
| 41 |
|
---|
| 42 |
|
---|
| 43 | typedef permutation_matrix<std::size_t> pmatrix;
|
---|
| 44 | // create a working copy of the input
|
---|
| 45 | matrix<RealType> A(m_cov);
|
---|
| 46 | // create a permutation matrix for the LU-factorization
|
---|
| 47 | pmatrix pm(A.size1());
|
---|
| 48 | // perform LU-factorization
|
---|
| 49 | int res = lu_factorize(A,pm);
|
---|
| 50 | if( res != 0 ) throw math_error("Pdf function : covariance matrix is a singular matrix");
|
---|
| 51 |
|
---|
| 52 | // create identity matrix of "inverse"
|
---|
| 53 | m_invcov = identity_matrix<double>(A.size1());
|
---|
| 54 | // backsubstitute to get the inverse
|
---|
| 55 | lu_substitute(A, pm, invcov);
|
---|
| 56 |
|
---|
| 57 | //compute determinant
|
---|
| 58 | m_detcov = 1.0;
|
---|
| 59 | for (std::size_t i=0; i < pm.size(); ++i) {
|
---|
| 60 | if (pm(i) != i)
|
---|
| 61 | m_detcov *= -1.0;
|
---|
| 62 | m_detcov *= A(i,i);
|
---|
| 63 | }
|
---|
| 64 |
|
---|
| 65 | }
|
---|
| 66 |
|
---|
| 67 | /*!
|
---|
| 68 | *\brief Get the inversion of covariance matrix
|
---|
| 69 | */
|
---|
| 70 | boost::numeric::ublas::matrix<RealType> invcov() const {
|
---|
| 71 | return m_invcov;
|
---|
| 72 | }
|
---|
| 73 |
|
---|
| 74 | /*!
|
---|
| 75 | *\brief Get determinant of covariance matrix
|
---|
| 76 | */
|
---|
| 77 | boost::numeric::ublas::matrix<RealType> detcov() const {
|
---|
| 78 | return m_detcov;
|
---|
| 79 | }
|
---|
| 80 |
|
---|
| 81 |
|
---|
| 82 | /*!
|
---|
| 83 | *\brief Get mean vector
|
---|
| 84 | */
|
---|
| 85 | boost::numeric::ublas::vector<RealType> mean() const {
|
---|
| 86 | return m_mean;
|
---|
| 87 | }
|
---|
| 88 |
|
---|
| 89 | /*!
|
---|
| 90 | *\brief Get covariance matrix
|
---|
| 91 | */
|
---|
| 92 | boost::numeric::ublas::matrix<RealType> cov() const {
|
---|
| 93 | return m_cov;
|
---|
| 94 | }
|
---|
| 95 |
|
---|
| 96 |
|
---|
| 97 |
|
---|
| 98 |
|
---|
| 99 | private :
|
---|
| 100 |
|
---|
| 101 |
|
---|
| 102 |
|
---|
| 103 |
|
---|
| 104 | /*!
|
---|
| 105 | *\brief mean vector
|
---|
| 106 | */
|
---|
| 107 | boost::numeric::ublas::vector<RealType> m_mean;
|
---|
| 108 |
|
---|
| 109 | /*!
|
---|
| 110 | *\brief covariance matrix
|
---|
| 111 | */
|
---|
| 112 | boost::numeric::ublas::matrix<RealType> m_cov;
|
---|
| 113 |
|
---|
| 114 |
|
---|
| 115 | /*!
|
---|
| 116 | *\brief covar
|
---|
| 117 | */
|
---|
| 118 | boost::numeric::ublas::matrix<RealType> m_invcov;
|
---|
| 119 |
|
---|
| 120 | /*!
|
---|
| 121 | *\brief cov matrix
|
---|
| 122 | */
|
---|
| 123 | double m_detcov;
|
---|
| 124 |
|
---|
| 125 | };
|
---|
| 126 |
|
---|
| 127 | typedef multivariate_normal_distribution<double> mvnormal;
|
---|
| 128 |
|
---|
| 129 |
|
---|
| 130 | /*!
|
---|
| 131 | *\fn inline RealType pdf(const multivariate_normal_distribution<RealType>& dist, const boost::numeric::ublas::vector<RealType> & x)
|
---|
| 132 | *\brief Compute probability density function for a multivariate normal distribution
|
---|
| 133 | */
|
---|
| 134 | template <class RealType>
|
---|
| 135 | inline RealType pdf(const multivariate_normal_distribution<RealType>& dist, const boost::numeric::ublas::vector<RealType> & x)
|
---|
| 136 | {
|
---|
| 137 |
|
---|
| 138 | boost::numeric::ublas::vector<RealType> mean = dist.mean();
|
---|
| 139 | boost::numeric::ublas::matrix<RealType> cov= dist.cov();
|
---|
| 140 | boost::numeric::ublas::matrix<RealType> invcov =dist.invcov();
|
---|
| 141 | double detcov= dist.detcov;
|
---|
| 142 |
|
---|
| 143 |
|
---|
| 144 | RealType exponent = - Dot(x-mean, invcov*(x-mean));
|
---|
| 145 | exponent /= 2 ;
|
---|
| 146 |
|
---|
| 147 | RealType result = std::exp(exponent);
|
---|
| 148 | result /= std::sqrt(pow(2*M_PI,mean.size())*std::abs(detcov));
|
---|
| 149 |
|
---|
| 150 | return result;
|
---|
| 151 | } // pdf
|
---|
| 152 |
|
---|
| 153 |
|
---|
| 154 | /*!
|
---|
| 155 | *\class multivariate_uniform_distribution
|
---|
| 156 | *\brief This clas describes a multivariate normal distribution
|
---|
| 157 | */
|
---|
| 158 | template < class RealType =double >
|
---|
| 159 | class multivariate_uniform_distribution{
|
---|
| 160 | public :
|
---|
| 161 |
|
---|
| 162 | /*!
|
---|
| 163 | *\brief Default constructor
|
---|
| 164 | */
|
---|
| 165 | multivariate_uniform_distribution(boost::numeric::ublas::vector<RealType> lower,boost::numeric::ublas::matrix<RealType> upper){
|
---|
| 166 |
|
---|
| 167 | if(upper.size()!= lower.size()) throw math_error("Multivariate uniform distribution : the upper vector and the loxer vector must have the same size");
|
---|
| 168 | for (size_t i=0;i<upper.size();i++)
|
---|
| 169 | if(lower[i]>upper[i]) throw math_error("Multivariate uniform distribution : the lower vector is not lower than upper vector");
|
---|
| 170 |
|
---|
| 171 | m_lower=lower;
|
---|
| 172 | m_upper=upper;
|
---|
| 173 | }
|
---|
| 174 | private :
|
---|
| 175 |
|
---|
| 176 | /*!
|
---|
| 177 | *\brief Lower vector
|
---|
| 178 | */
|
---|
| 179 | boost::numeric::ublas::vector<RealType> m_lower;
|
---|
| 180 |
|
---|
| 181 | /*!
|
---|
| 182 | *\brief Upper vector
|
---|
| 183 | */
|
---|
| 184 | boost::numeric::ublas::matrix<RealType> m_upper;
|
---|
| 185 |
|
---|
| 186 | };
|
---|
| 187 |
|
---|
| 188 | typedef multivariate_normal_distribution<double> mvuniform;
|
---|
| 189 |
|
---|
| 190 |
|
---|
| 191 | /*!
|
---|
| 192 | *\fn inline RealType pdf(const multivariate_uniform_distribution<RealType>& dist, const boost::numeric::ublas::vector<RealType> & x)
|
---|
| 193 | *\brief Compute probability density function for a multivariate uniform distribution
|
---|
| 194 | */
|
---|
| 195 | template <class RealType>
|
---|
| 196 | inline RealType pdf(const multivariate_uniform_distribution<RealType>& dist, const boost::numeric::ublas::vector<RealType> & x)
|
---|
| 197 | {
|
---|
| 198 |
|
---|
| 199 | double result =1;
|
---|
| 200 |
|
---|
| 201 | boost::numeric::ublas::vector<RealType> upper = dist.upper();
|
---|
| 202 | boost::numeric::ublas::matrix<RealType> lower= dist.lower();
|
---|
| 203 |
|
---|
| 204 |
|
---|
| 205 | for(size_t i=0;i<upper.size();i++){
|
---|
| 206 | if(x<upper[i] && x[i]>lower[i]) result/=(upper[i]-lower[i]);
|
---|
| 207 | else {result=0; break;}
|
---|
| 208 | }
|
---|
| 209 |
|
---|
| 210 | return result;
|
---|
| 211 | } // pdf
|
---|
| 212 |
|
---|
| 213 | };
|
---|
| 214 | };
|
---|
| 215 | #endif
|
---|