gaussian.hpp

Go to the documentation of this file.

#ifndef SSMPACK_DISTRIBUTION_GAUSSIAN_HPP
#define SSMPACK_DISTRIBUTION_GAUSSIAN_HPP

#include "ssmkit/random/generator.hpp"

#include <armadillo>

#include <cmath>

namespace ssmkit {
namespace distribution {

class Gaussian {

 public:
  using TParameterVAR = std::tuple<arma::vec, arma::mat>;

 private:
  arma::vec mean_;
  arma::mat covariance_;
  std::normal_distribution<double> normal_;
  static constexpr double pi = 3.1415926535897;
  arma::mat inv_cov_;
  double part_;
  arma::mat chol_dec_;
  int dim_;

 private:
  void calcDistConstants() {
    dim_ = mean_.n_rows;
    inv_cov_ = arma::inv(covariance_);
    // calculate partition function
    double det_cov = arma::det(covariance_);
    double den_pi = 1 / std::sqrt(std::pow(2 * pi, dim_));
    part_ = den_pi * (1 / std::sqrt(det_cov));
    // Cholesky decomposition
    chol_dec_ = arma::chol(covariance_, "lower");
  }

 public:
  Gaussian() = delete;
  Gaussian(int dim) : Gaussian(arma::zeros(dim), arma::eye(dim, dim)) {}
  
  Gaussian(arma::vec mean, arma::mat covariance)
      : mean_(std::move(mean)), covariance_(std::move(covariance)) {
    calcDistConstants();
  }

  arma::vec random() {
    arma::vec rnd(dim_); // how much overload the vector construction has?
    rnd.imbue(
        [&]() { return normal_(random::Generator::get().getGenerator()); });
    return mean_ + chol_dec_ * rnd;
  }

  double likelihood(const arma::vec &rv) const {
    const auto diff = rv - mean_;
    const arma::vec tmp = diff.t() * inv_cov_ * diff;
    const double expt = -tmp(0) / 2;
    return part_ * std::exp(expt);
  }

  Gaussian &parameterize(const TParameterVAR &parameters) {
    return parameterize(std::get<0>(parameters), std::get<1>(parameters));
  }

  Gaussian &parameterize(const arma::vec &mean,
                         const arma::mat &covariance) {
    mean_ = mean;
    covariance_ = covariance;
    calcDistConstants();
    return (*this);
  }
  const arma::vec& getMean() const { return mean_; }
  const arma::mat& getCovariance() const { return covariance_; }
};

} // namespace distribution
} // namespace ssmkit

#endif // SSMPACK_DISTRIBUTION_GAUSSIAN_HPP