MultiVariateNormalDistribution.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Implements a multi-variate normal distribution with zero mean.
5  *
6  *
7  *
8  * \author T.Voss
9  * \date 2010
10  *
11  *
12  * \par Copyright 1995-2015 Shark Development Team
13  *
14  * <BR><HR>
15  * This file is part of Shark.
16  * <http://image.diku.dk/shark/>
17  *
18  * Shark is free software: you can redistribute it and/or modify
19  * it under the terms of the GNU Lesser General Public License as published
20  * by the Free Software Foundation, either version 3 of the License, or
21  * (at your option) any later version.
22  *
23  * Shark is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26  * GNU Lesser General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License
29  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
30  *
31  */
32 #ifndef SHARK_STATISTICS_MULTIVARIATENORMALDISTRIBUTION_H
33 #define SHARK_STATISTICS_MULTIVARIATENORMALDISTRIBUTION_H
34 
36 #include <shark/LinAlg/Cholesky.h>
37 #include <shark/Rng/GlobalRng.h>
38 
39 namespace shark {
40 
41 /// \brief Implements a multi-variate normal distribution with zero mean.
43 public:
44  ///\brief Result type of a sampling operation.
45  ///
46  /// The first element is the result of sampling this distribution, the
47  /// second element is the original standard-normally distributed vector drawn
48  /// for sampling purposes.
49  typedef std::pair<RealVector,RealVector> result_type;
50 
51  /// \brief Constructor
52  /// \param [in] Sigma covariance matrix
53  MultiVariateNormalDistribution( RealMatrix const& Sigma ) {
54  m_covarianceMatrix = Sigma;
55  update();
56  }
57 
59 
60  /// \brief Stores/Restores the distribution from the supplied archive.
61  /// \param [in,out] ar The archive to read from/write to.
62  /// \param [in] version Currently unused.
63  template<typename Archive>
64  void serialize( Archive & ar, const std::size_t version ) {
65  ar & BOOST_SERIALIZATION_NVP( m_covarianceMatrix );
66  ar & BOOST_SERIALIZATION_NVP( m_eigenVectors );
67  ar & BOOST_SERIALIZATION_NVP( m_eigenValues );
68  }
69 
70  /// \brief Resizes the distribution. Updates both eigenvectors and eigenvalues.
71  /// \param [in] size The new size of the distribution
72  void resize( std::size_t size ) {
73  m_covarianceMatrix = blas::identity_matrix<double>( size );
74  m_eigenValues = blas::repeat(1.0,size);
75  m_eigenVectors = blas::identity_matrix<double>( size );
76  }
77 
78  /// \brief Accesses the covariance matrix defining the distribution.
79  RealMatrix const& covarianceMatrix() const {
80  return m_covarianceMatrix;
81  }
82 
83  /// \brief Accesses a mutable reference to the covariance matrix
84  /// defining the distribution. Allows for l-value semantics.
85  ///
86  /// ATTENTION: If the reference is altered, update needs to be called manually.
87  RealMatrix& covarianceMatrix() {
88  return m_covarianceMatrix;
89  }
90 
91  /// \brief Sets the covariance matrix and updates the internal variables. This is expensive
92  void setCovarianceMatrix(RealMatrix const& matrix){
93  covarianceMatrix() = matrix;
94  update();
95  }
96 
97  /// \brief Accesses an immutable reference to the eigenvectors of the covariance matrix.
98  const RealMatrix & eigenVectors() const {
99  return m_eigenVectors;
100  }
101 
102  /// \brief Accesses an immutable reference to the eigenvalues of the covariance matrix.
103  RealVector const& eigenValues() const {
104  return m_eigenValues;
105  }
106 
107  /// \brief Accesses a reference to the eigenvalues of the covariance matrix.
108  RealVector& eigenValues(){
109  return m_eigenValues;
110  }
111 
112  /// \brief Samples the distribution.
113  result_type operator()() const {
114  RealVector result( m_eigenValues.size(), 0. );
115  RealVector z( m_eigenValues.size() );
116 
117  for( std::size_t i = 0; i < result.size(); i++ ) {
118  z( i ) = Rng::gauss( 0., 1. );
119  }
120 
121  for( std::size_t i = 0; i < result.size(); i++ )
122  for( std::size_t j = 0; j < result.size(); j++ )
123  result( i ) += m_eigenVectors( i, j ) * std::sqrt( std::abs( m_eigenValues(j) ) ) * z( j );
124 
125  return( std::make_pair( result, z ) );
126  }
127 
128  /// \brief Calculates the evd of the current covariance matrix.
129  void update() {
130  eigensymm( m_covarianceMatrix, m_eigenVectors, m_eigenValues );
131  }
132 
133 private:
134  RealMatrix m_covarianceMatrix; ///< Covariance matrix of the mutation distribution.
135  RealMatrix m_eigenVectors; ///< Eigenvectors of the covariance matrix.
136  RealVector m_eigenValues; ///< Eigenvalues of the covariance matrix.
137 };
138 
139 /// \brief Multivariate normal distribution with zero mean using a cholesky decomposition
141 public:
142  /// \brief Result type of a sampling operation.
143  ///
144  /// The first element is the result of sampling this distribution, the
145  /// second element is the original standard-normally distributed vector drawn
146  /// for sampling purposes.
147  typedef std::pair<RealVector,RealVector> result_type;
148 
149  /// \brief Constructor
150  /// \param [in] covariance covariance matrix
151  /// \param triangular Is the choleksy factor triangular?
152  MultiVariateNormalDistributionCholesky( RealMatrix const& covariance, bool triangular=false )
153  :m_triangular(false){
154  setCovarianceMatrix(covariance);
155  }
156 
157  MultiVariateNormalDistributionCholesky(bool triangular=false):m_triangular(triangular){}
158 
159  /// \brief Stores/Restores the distribution from the supplied archive.
160  ///\param [in,out] ar The archive to read from/write to.
161  ///\param [in] version Currently unused.
162  template<typename Archive>
163  void serialize( Archive & ar, const std::size_t version ) {
164  ar & BOOST_SERIALIZATION_NVP( m_lowerCholesky);
165  }
166 
167  /// \brief Resizes the distribution. Updates both eigenvectors and eigenvalues.
168  /// \param [in] size The new size of the distribution
169  void resize( std::size_t size ) {
170  m_lowerCholesky = blas::identity_matrix<double>( size );
171  }
172 
173  /// \brief Returns the size of the created vectors
174  std::size_t size()const{
175  return m_lowerCholesky.size1();
176  }
177 
178  /// \brief Sets the new covariance matrix by computing the new cholesky dcomposition
179  void setCovarianceMatrix(RealMatrix const& matrix){
180  choleskyDecomposition(matrix,m_lowerCholesky);
181  }
182 
183  /// \brief Returns the lower cholsky factor.
185  return m_lowerCholesky;
186  }
187 
188  /// \brief Returns the lower cholesky factor.
190  return m_lowerCholesky;
191  }
192 
193  void rankOneUpdate(double alpha, double beta, RealVector const& v){
194  choleskyUpdate(m_lowerCholesky,v,alpha,beta);
195  }
196 
197 
198  template<class Vector1, class Vector2>
199  void generate(Vector1& y, Vector2& z)const{
200  z.resize(size());
201  y.resize(size());
202 
203  for( std::size_t i = 0; i != size(); i++ ) {
204  z( i ) = Rng::gauss( 0, 1 );
205  }
206 
207  if(m_triangular && size() > 400){
208  y=z;
209  blas::triangular_prod<blas::lower>(m_lowerCholesky,y);
210  }else{
211  noalias(y) = prod(m_lowerCholesky,z);
212  }
213  }
214 
215  /// \brief Samples the distribution.
216  ///
217  /// Returns a vector pair (y,z) where y=Lz and, L is the lower cholesky factor and z is a vector
218  /// of normally distributed numbers. Thus y is the real sampled point.
219  result_type operator()() const {
220  result_type result;
221  RealVector& z = result.second;
222  RealVector& y = result.first;
223  generate(y,z);
224  return result;
225  }
226 
227 private:
228  blas::matrix<double,blas::column_major> m_lowerCholesky; ///< The lower cholesky factor (actually any is okay as long as it is the left)
229  bool m_triangular;
230 };
231 
232 
233 }
234 
235 #endif