forked from mlpack/mlpack
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gaussian_distribution.hpp
188 lines (156 loc) · 5.26 KB
/
gaussian_distribution.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
/**
* @file gaussian_distribution.hpp
* @author Ryan Curtin
* @author Michael Fox
*
* Implementation of the Gaussian distribution.
*/
#ifndef __MLPACK_CORE_DISTRIBUTIONS_GAUSSIAN_DISTRIBUTION_HPP
#define __MLPACK_CORE_DISTRIBUTIONS_GAUSSIAN_DISTRIBUTION_HPP
#include <mlpack/core.hpp>
namespace mlpack {
namespace distribution {
/**
* A single multivariate Gaussian distribution.
*/
class GaussianDistribution
{
private:
//! Mean of the distribution.
arma::vec mean;
//! Positive definite covariance of the distribution.
arma::mat covariance;
//! Lower triangular factor of cov (e.g. cov = LL^T).
arma::mat covLower;
//! Cached inverse of covariance.
arma::mat invCov;
//! Cached logdet(cov).
double logDetCov;
//! log(2pi)
static const constexpr double log2pi = 1.83787706640934533908193770912475883;
public:
/**
* Default constructor, which creates a Gaussian with zero dimension.
*/
GaussianDistribution() { /* nothing to do */ }
/**
* Create a Gaussian distribution with zero mean and identity covariance with
* the given dimensionality.
*/
GaussianDistribution(const size_t dimension) :
mean(arma::zeros<arma::vec>(dimension)),
covariance(arma::eye<arma::mat>(dimension, dimension)),
covLower(arma::eye<arma::mat>(dimension, dimension)),
invCov(arma::eye<arma::mat>(dimension, dimension)),
logDetCov(0)
{ /* Nothing to do. */ }
/**
* Create a Gaussian distribution with the given mean and covariance.
*
* covariance is expected to be positive definite.
*/
GaussianDistribution(const arma::vec& mean, const arma::mat& covariance);
// TODO(stephentu): do we want a (arma::vec&&, arma::mat&&) ctor?
//! Return the dimensionality of this distribution.
size_t Dimensionality() const { return mean.n_elem; }
/**
* Return the probability of the given observation.
*/
double Probability(const arma::vec& observation) const
{
return exp(LogProbability(observation));
}
/**
* Return the log probability of the given observation.
*/
double LogProbability(const arma::vec& observation) const;
/**
* Calculates the multivariate Gaussian probability density function for each
* data point (column) in the given matrix
*
* @param x List of observations.
* @param probabilities Output probabilities for each input observation.
*/
void Probability(const arma::mat& x, arma::vec& probabilities) const
{
arma::vec logProbabilities;
LogProbability(x, logProbabilities);
probabilities = arma::exp(logProbabilities);
}
void LogProbability(const arma::mat& x, arma::vec& logProbabilities) const;
/**
* Return a randomly generated observation according to the probability
* distribution defined by this object.
*
* @return Random observation from this Gaussian distribution.
*/
arma::vec Random() const;
/**
* Estimate the Gaussian distribution directly from the given observations.
*
* @param observations List of observations.
*/
void Estimate(const arma::mat& observations);
/**
* Estimate the Gaussian distribution from the given observations, taking into
* account the probability of each observation actually being from this
* distribution.
*/
void Estimate(const arma::mat& observations,
const arma::vec& probabilities);
/**
* Return the mean.
*/
const arma::vec& Mean() const { return mean; }
/**
* Return a modifiable copy of the mean.
*/
arma::vec& Mean() { return mean; }
/**
* Return the covariance matrix.
*/
const arma::mat& Covariance() const { return covariance; }
/**
* Set the covariance.
*/
void Covariance(const arma::mat& covariance);
void Covariance(arma::mat&& covariance);
/**
* Returns a string representation of this object.
*/
std::string ToString() const;
/*
* Save to or Load from SaveRestoreUtility
*/
void Save(util::SaveRestoreUtility& n) const;
void Load(const util::SaveRestoreUtility& n);
static std::string const Type() { return "GaussianDistribution"; }
private:
void FactorCovariance();
};
/**
* Calculates the multivariate Gaussian log probability density function for each
* data point (column) in the given matrix
*
* @param x List of observations.
* @param probabilities Output log probabilities for each input observation.
*/
inline void GaussianDistribution::LogProbability(const arma::mat& x,
arma::vec& logProbabilities) const
{
// Column i of 'diffs' is the difference between x.col(i) and the mean.
arma::mat diffs = x - (mean * arma::ones<arma::rowvec>(x.n_cols));
// Now, we only want to calculate the diagonal elements of (diffs' * cov^-1 *
// diffs). We just don't need any of the other elements. We can calculate
// the right hand part of the equation (instead of the left side) so that
// later we are referencing columns, not rows -- that is faster.
const arma::mat rhs = -0.5 * invCov * diffs;
arma::vec logExponents(diffs.n_cols); // We will now fill this.
for (size_t i = 0; i < diffs.n_cols; i++)
logExponents(i) = accu(diffs.unsafe_col(i) % rhs.unsafe_col(i));
const size_t k = x.n_rows;
logProbabilities = -0.5 * k * log2pi - 0.5 * logDetCov + logExponents;
}
}; // namespace distribution
}; // namespace mlpack
#endif