rtmvnorm.h

Go to the documentation of this file.

/* 
 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
 * and Daniel Pemstein.  All Rights Reserved.
 *
 * This program is free software; you can redistribute it and/or
 * modify under the terms of the GNU General Public License as
 * published by Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.  See the text files
 * COPYING and LICENSE, distributed with this source code, for further
 * information.
 * --------------------------------------------------------------------
 *  scythestat/rng/rtmvnorm.h
 *
 */

#ifndef SCYTHE_RTMVNORM_H
#define SCYTHE_RTMVNORM_H

#include <iostream>
#include <cmath>

#ifdef SCYTHE_COMPILE_DIRECT
#include "matrix.h"
#include "rng.h"
#include "error.h"
#include "algorithm.h"
#include "ide.h"
#else
#include "scythestat/matrix.h"
#include "scythestat/rng.h"
#include "scythestat/error.h"
#include "scythestat/algorithm.h"
#include "scythestat/ide.h"
#endif
namespace scythe
{
  /* Truncated Multivariate Normal Distribution by Gibbs sampling
   * (Geweke 1991).  This is a functor that allows one to
   * initialize---and optionally burn in---a sampler for a given
   * truncated multivariate normal distribution on construction
   * and then make (optionally thinned) draws with calls to the ()
   * operator.
   *
   */
  template <class RNGTYPE>
  class rtmvnorm {
    public:

      template <matrix_order PO1, matrix_style PS1, matrix_order PO2,
                matrix_style PS2, matrix_order PO3, matrix_style PS3,
                matrix_order PO4, matrix_style PS5, matrix_order PO5,
                matrix_style PS4>
      rtmvnorm (const Matrix<double, PO1,PS1>& mu,
                const Matrix<double, PO2, PS2>& sigma,
                const Matrix<double, PO3, PS3>& D,
                const Matrix<double, PO4, PS4>& a,
                const Matrix<double, PO5, PS5>& b, rng<RNGTYPE>& generator, 
                unsigned int burnin = 0, unsigned int thin = 1,
                bool preinvertedD = false) 
      : mu_ (mu), C_ (mu.rows(), mu.rows(), false), 
        h_ (mu.rows(), 1, false), z_ (mu.rows(), 1, true, 0), 
        generator_ (generator), n_ (mu.rows()), thin_ (thin), iter_ (0)
      {
        SCYTHE_CHECK_10(thin == 0, scythe_invalid_arg,
            "thin must be >= 1");
        SCYTHE_CHECK_10(! mu.isColVector(), scythe_dimension_error,
            "mu not column vector");
        SCYTHE_CHECK_10(!  sigma.isSquare(), scythe_dimension_error,
            "sigma not square");
        SCYTHE_CHECK_10(!  D.isSquare(), scythe_dimension_error,
            "D not square");
        SCYTHE_CHECK_10(!  a.isColVector(), scythe_dimension_error,
            "a not column vector");
        SCYTHE_CHECK_10(!  b.isColVector(), scythe_dimension_error,
            "b not column vector");
        SCYTHE_CHECK_10(sigma.rows() != n_ || D.rows() != n_ || 
            a.rows() != n_ || b.rows() != n_, scythe_conformation_error,
            "mu, sigma, D, a, and b not conformable");

        // TODO will D * sigma * t(D) always be positive definite,
        // allowing us to use the faster invpd?
        if (preinvertedD)
          Dinv_ = D;
        else
          Dinv_ = inv(D);
        Matrix<> Tinv = inv(D * sigma * t(D));
        alpha_ = a - D * mu;
        beta_ =  b - D * mu;

        // Check truncation bounds
        if (SCYTHE_DEBUG > 0) {
          for (unsigned int i = 0; i < n_; ++i) {
            SCYTHE_CHECK(alpha_(i) >= beta_(i), scythe_invalid_arg,
                "Truncation bound " << i 
                << " not logically consistent");
          }
        }

        // Precompute some stuff (see Geweke 1991 pg 7).
        for (unsigned int i = 0; i < n_; ++i) {
          C_(i, _) = -(1 / Tinv(i, i)) % Tinv(i, _);
          C_(i, i) = 0; // not really clever but probably too clever
          h_(i) = std::sqrt(1 / Tinv(i, i));
          SCYTHE_CHECK_30(std::isnan(h_(i)), scythe_invalid_arg,
              "sigma is not positive definite");
        }

        // Do burnin
        for (unsigned int i = 0; i < burnin; ++i)
          sample ();
      }

      template <matrix_order O, matrix_style S>
      Matrix<double, O, S> operator() ()
      {
        do { sample (); } while (iter_ % thin_ != 0);

        return (mu_ + Dinv_ * z_);
      }

      Matrix<double,Col,Concrete> operator() ()
      {
        return operator()<Col, Concrete>();
      }

    protected:
      /* Does one step of the Gibbs sampler (see Geweke 1991 p 6) */
      void sample ()
      {
        double czsum;
        double above;
        double below;
        for (unsigned int i = 0; i < n_; ++i) {

          // Calculate sum_{j \ne i} c_{ij} z_{j}
          czsum = 0;
          for (unsigned int j = 0; j < n_; ++j) {
            if (i == j) continue;
            czsum += C_(i, j) * z_(j);
          }

          // Calc truncation of conditional univariate std normal
          below = (alpha_(i) - czsum) / h_(i);
          above = (beta_(i) - czsum) / h_(i);
          
          // Draw random variate z_i
          z_(i) = h_(i);
          if (above == std::numeric_limits<double>::infinity()){
            if (below == -std::numeric_limits<double>::infinity())
              z_(i) *= generator_.rnorm(0, 1); // untruncated
            else
              z_(i) *= generator_.rtbnorm_combo(0, 1, below);
          } else if (below == 
              -std::numeric_limits<double>::infinity())
            z_(i) *= generator_.rtanorm_combo(0, 1, above);
          else
            z_(i) *= generator_.rtnorm_combo(0, 1, below, above);

          z_(i) += czsum;
        }

        ++iter_;
      }

      /* Instance variables */
      // Various reused computation matrices with names from
      // Geweke 1991.
      Matrix<> mu_;   Matrix<> Dinv_; 
      Matrix<> C_; Matrix<> alpha_; Matrix<> beta_; Matrix<> h_; 

      Matrix<> z_; // The current draw of the posterior

      rng<RNGTYPE>& generator_; // Refernce to random number generator

      unsigned int n_;  // The dimension of the distribution
      unsigned int thin_; // thinning parameter
      unsigned int iter_; // The current post-burnin iteration
  };
} // end namespace scythe
#endif