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