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() : logDetCov(0.0) { /* nothing to do */ }

  GaussianDistribution(const size_t dimension) :
      mean(arma::zeros(dimension)),
      covariance(arma::eye(dimension, dimension)),
      covLower(arma::eye(dimension, dimension)),
      invCov(arma::eye(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
  {
    probabilities.set_size(x.n_cols);
    for (size_t i = 0; i < x.n_cols; ++i)
    {
      probabilities(i) = Probability(x.unsafe_col(i));
    }
  }

  void 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;
    diffs.each_col() -= mean;
    // 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));

    logProbabilities = -0.5 * x.n_rows * log2pi - 0.5 * logDetCov +
      logExponents;
  }

  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);

  const arma::mat& InvCov() const { return invCov; }

  double LogDetCov() const { return logDetCov; }

  template<typename Archive>
  void serialize(Archive& ar, const unsigned int /* version */)
  {
    // We just need to serialize each of the members.
    ar & BOOST_SERIALIZATION_NVP(mean);
    ar & BOOST_SERIALIZATION_NVP(covariance);
    ar & BOOST_SERIALIZATION_NVP(covLower);
    ar & BOOST_SERIALIZATION_NVP(invCov);
    ar & BOOST_SERIALIZATION_NVP(logDetCov);
  }

 private:
  void FactorCovariance();
};

} // namespace distribution
} // namespace mlpack

#endif