source: pacpusframework/trunk/include/Pacpus/PacpusTools/math/distributions.hpp

Last change on this file was 89, checked in by morasjul, 11 years ago

PACPUS 2.0 Beta deployed in new branch

Major changes:
-Add communication interface between components
-Add examples for communications interface (TestComponents)
-Move to Qt5 support

  • Property svn:executable set to *
File size: 5.8 KB
Line 
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
20namespace math{
21namespace distributions{
22
23
24template <class RealType>
25inline RealType pdf(const boost::math::normal_distribution<RealType> & dist , boost::numeric::ublas::vector<RealType> & v){
26 return pdf(dist,v[0]);
27}
28
29template <class RealType>
30inline 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 */
39template < class RealType =double >
40class multivariate_normal_distribution{
41public :
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
112private :
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
140typedef 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 */
147template <class RealType>
148inline 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 */
171template < class RealType =double >
172class multivariate_uniform_distribution{
173public :
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
201typedef 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 */
208template <class RealType>
209inline 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
Note: See TracBrowser for help on using the repository browser.