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