gaussian_distribution.hpp

Go to the documentation of this file.

#ifndef MLPACK_CORE_DISTRIBUTIONS_GAUSSIAN_DISTRIBUTION_HPP
#define MLPACK_CORE_DISTRIBUTIONS_GAUSSIAN_DISTRIBUTION_HPP

#include <mlpack/prereqs.hpp>

namespace mlpack {
namespace distribution {

class GaussianDistribution
{
 private:
  arma::vec mean;
  arma::mat covariance;
  arma::mat covLower;
  arma::mat invCov;
  double logDetCov;

  static const constexpr double log2pi = 1.83787706640934533908193770912475883;

 public:
  GaussianDistribution() { /* nothing to do */ }

  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. */ }

  GaussianDistribution(const arma::vec& mean, const arma::mat& covariance);

  // TODO(stephentu): do we want a (arma::vec&&, arma::mat&&) ctor?

  size_t Dimensionality() const { return mean.n_elem; }

  double Probability(const arma::vec& observation) const
  {
    return exp(LogProbability(observation));
  }

  double LogProbability(const arma::vec& observation) const;

  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;

  arma::vec Random() const;

  void Train(const arma::mat& observations);

  void Train(const arma::mat& observations,
             const arma::vec& probabilities);

  const arma::vec& Mean() const { return mean; }

  arma::vec& Mean() { return mean; }

  const arma::mat& Covariance() const { return covariance; }

  void Covariance(const arma::mat& covariance);

  void Covariance(arma::mat&& covariance);

  template<typename Archive>
  void Serialize(Archive& ar, const unsigned int /* version */)
  {
    using data::CreateNVP;

    // We just need to serialize each of the members.
    ar & CreateNVP(mean, "mean");
    ar & CreateNVP(covariance, "covariance");
    ar & CreateNVP(covLower, "covLower");
    ar & CreateNVP(invCov, "invCov");
    ar & CreateNVP(logDetCov, "logDetCov");
  }

 private:
  void FactorCovariance();
};

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