Stan  2.10.0
probability, sampling & optimization
normal_meanfield.hpp
Go to the documentation of this file.
1 #ifndef STAN_VARIATIONAL_NORMAL_MEANFIELD_HPP
2 #define STAN_VARIATIONAL_NORMAL_MEANFIELD_HPP
3 
5 #include <stan/math/prim/mat/fun/Eigen.hpp>
6 #include <stan/math/prim/scal/fun/constants.hpp>
7 #include <stan/math/prim/mat/meta/get.hpp>
8 #include <stan/math/prim/arr/meta/get.hpp>
9 #include <stan/math/prim/mat/meta/length.hpp>
10 #include <stan/math/prim/mat/meta/is_vector.hpp>
11 #include <stan/math/prim/mat/meta/is_vector_like.hpp>
12 #include <stan/math/prim/mat/fun/value_of_rec.hpp>
13 #include <stan/math/prim/scal/prob/normal_rng.hpp>
14 #include <stan/math/prim/scal/err/check_finite.hpp>
15 #include <stan/math/prim/scal/err/check_not_nan.hpp>
16 #include <stan/math/prim/scal/err/check_positive.hpp>
17 #include <stan/math/prim/scal/err/check_size_match.hpp>
18 #include <stan/math/prim/scal/err/domain_error.hpp>
19 #include <stan/model/util.hpp>
21 #include <algorithm>
22 #include <ostream>
23 #include <vector>
24 
25 namespace stan {
26 
27  namespace variational {
28 
33  class normal_meanfield : public base_family {
34  private:
38  Eigen::VectorXd mu_;
39 
43  Eigen::VectorXd omega_;
44 
48  const int dimension_;
49 
50  public:
58  explicit normal_meanfield(size_t dimension)
59  : mu_(Eigen::VectorXd::Zero(dimension)),
60  omega_(Eigen::VectorXd::Zero(dimension)),
61  dimension_(dimension) {
62  }
63 
71  explicit normal_meanfield(const Eigen::VectorXd& cont_params)
72  : mu_(cont_params),
73  omega_(Eigen::VectorXd::Zero(cont_params.size())),
74  dimension_(cont_params.size()) {
75  }
76 
87  normal_meanfield(const Eigen::VectorXd& mu,
88  const Eigen::VectorXd& omega)
89  : mu_(mu), omega_(omega), dimension_(mu.size()) {
90  static const char* function =
91  "stan::variational::normal_meanfield";
92  stan::math::check_size_match(function,
93  "Dimension of mean vector", dimension_,
94  "Dimension of log std vector", omega_.size() );
95  stan::math::check_not_nan(function, "Mean vector", mu_);
96  stan::math::check_not_nan(function, "Log std vector", omega_);
97  }
98 
102  int dimension() const { return dimension_; }
103 
107  const Eigen::VectorXd& mu() const { return mu_; }
108 
112  const Eigen::VectorXd& omega() const { return omega_; }
113 
122  void set_mu(const Eigen::VectorXd& mu) {
123  static const char* function =
124  "stan::variational::normal_meanfield::set_mu";
125 
126  stan::math::check_size_match(function,
127  "Dimension of input vector", mu.size(),
128  "Dimension of current vector", dimension_);
129  stan::math::check_not_nan(function, "Input vector", mu);
130  mu_ = mu;
131  }
132 
142  void set_omega(const Eigen::VectorXd& omega) {
143  static const char* function =
144  "stan::variational::normal_meanfield::set_omega";
145 
146  stan::math::check_size_match(function,
147  "Dimension of input vector", omega.size(),
148  "Dimension of current vector", dimension_);
149  stan::math::check_not_nan(function, "Input vector", omega);
150  omega_ = omega;
151  }
152 
157  void set_to_zero() {
158  mu_ = Eigen::VectorXd::Zero(dimension_);
159  omega_ = Eigen::VectorXd::Zero(dimension_);
160  }
161 
169  return normal_meanfield(Eigen::VectorXd(mu_.array().square()),
170  Eigen::VectorXd(omega_.array().square()));
171  }
172 
184  return normal_meanfield(Eigen::VectorXd(mu_.array().sqrt()),
185  Eigen::VectorXd(omega_.array().sqrt()));
186  }
187 
200  static const char* function =
201  "stan::variational::normal_meanfield::operator=";
202  stan::math::check_size_match(function,
203  "Dimension of lhs", dimension_,
204  "Dimension of rhs", rhs.dimension());
205  mu_ = rhs.mu();
206  omega_ = rhs.omega();
207  return *this;
208  }
209 
222  static const char* function =
223  "stan::variational::normal_meanfield::operator+=";
224  stan::math::check_size_match(function,
225  "Dimension of lhs", dimension_,
226  "Dimension of rhs", rhs.dimension());
227  mu_ += rhs.mu();
228  omega_ += rhs.omega();
229  return *this;
230  }
231 
245  static const char* function =
246  "stan::variational::normal_meanfield::operator/=";
247  stan::math::check_size_match(function,
248  "Dimension of lhs", dimension_,
249  "Dimension of rhs", rhs.dimension());
250  mu_.array() /= rhs.mu().array();
251  omega_.array() /= rhs.omega().array();
252  return *this;
253  }
254 
266  normal_meanfield& operator+=(double scalar) {
267  mu_.array() += scalar;
268  omega_.array() += scalar;
269  return *this;
270  }
271 
284  normal_meanfield& operator*=(double scalar) {
285  mu_ *= scalar;
286  omega_ *= scalar;
287  return *this;
288  }
289 
297  const Eigen::VectorXd& mean() const {
298  return mu();
299  }
300 
311  double entropy() const {
312  return 0.5 * static_cast<double>(dimension_) *
313  (1.0 + stan::math::LOG_TWO_PI) + omega_.sum();
314  }
315 
328  Eigen::VectorXd transform(const Eigen::VectorXd& eta) const {
329  static const char* function =
330  "stan::variational::normal_meanfield::transform";
331  stan::math::check_size_match(function,
332  "Dimension of mean vector", dimension_,
333  "Dimension of input vector", eta.size() );
334  stan::math::check_not_nan(function, "Input vector", eta);
335  // exp(omega) * eta + mu
336  return eta.array().cwiseProduct(omega_.array().exp()) + mu_.array();
337  }
338 
349  template <class BaseRNG>
350  void sample(BaseRNG& rng, Eigen::VectorXd& eta) const {
351  // Draw from standard normal and transform to real-coordinate space
352  for (int d = 0; d < dimension_; ++d)
353  eta(d) = stan::math::normal_rng(0, 1, rng);
354  eta = transform(eta);
355  }
356 
375  template <class M, class BaseRNG>
376  void calc_grad(normal_meanfield& elbo_grad,
377  M& m,
378  Eigen::VectorXd& cont_params,
379  int n_monte_carlo_grad,
380  BaseRNG& rng,
382  const {
383  static const char* function =
384  "stan::variational::normal_meanfield::calc_grad";
385 
386  stan::math::check_size_match(function,
387  "Dimension of elbo_grad", elbo_grad.dimension(),
388  "Dimension of variational q", dimension_);
389  stan::math::check_size_match(function,
390  "Dimension of variational q", dimension_,
391  "Dimension of variables in model", cont_params.size());
392 
393  Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
394  Eigen::VectorXd omega_grad = Eigen::VectorXd::Zero(dimension_);
395  double tmp_lp = 0.0;
396  Eigen::VectorXd tmp_mu_grad = Eigen::VectorXd::Zero(dimension_);
397  Eigen::VectorXd eta = Eigen::VectorXd::Zero(dimension_);
398  Eigen::VectorXd zeta = Eigen::VectorXd::Zero(dimension_);
399 
400  // Naive Monte Carlo integration
401  static const int n_retries = 10;
402  for (int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
403  // Draw from standard normal and transform to real-coordinate space
404  for (int d = 0; d < dimension_; ++d)
405  eta(d) = stan::math::normal_rng(0, 1, rng);
406  zeta = transform(eta);
407  try {
408  std::stringstream ss;
409  stan::model::gradient(m, zeta, tmp_lp, tmp_mu_grad, &ss);
410  if (ss.str().length() > 0)
411  message_writer(ss.str());
412  stan::math::check_finite(function, "Gradient of mu", tmp_mu_grad);
413  mu_grad += tmp_mu_grad;
414  omega_grad.array() += tmp_mu_grad.array().cwiseProduct(eta.array());
415  ++i;
416  } catch (const std::exception& e) {
417  ++n_monte_carlo_drop;
418  if (n_monte_carlo_drop >= n_retries * n_monte_carlo_grad) {
419  const char* name = "The number of dropped evaluations";
420  const char* msg1 = "has reached its maximum amount (";
421  int y = n_retries * n_monte_carlo_grad;
422  const char* msg2 = "). Your model may be either severely "
423  "ill-conditioned or misspecified.";
424  stan::math::domain_error(function, name, y, msg1, msg2);
425  }
426  }
427  }
428  mu_grad /= static_cast<double>(n_monte_carlo_grad);
429  omega_grad /= static_cast<double>(n_monte_carlo_grad);
430 
431  omega_grad.array()
432  = omega_grad.array().cwiseProduct(omega_.array().exp());
433 
434  omega_grad.array() += 1.0; // add entropy gradient (unit)
435 
436  elbo_grad.set_mu(mu_grad);
437  elbo_grad.set_omega(omega_grad);
438  }
439  };
440 
451  const normal_meanfield& rhs) {
452  return lhs += rhs;
453  }
454 
465  const normal_meanfield& rhs) {
466  return lhs /= rhs;
467  }
468 
479  return rhs += scalar;
480  }
481 
492  return rhs *= scalar;
493  }
494 
495  }
496 }
497 #endif
void set_omega(const Eigen::VectorXd &omega)
Set the log standard deviation vector to the specified value.
const Eigen::VectorXd & mean() const
Returns the mean vector for this approximation.
base_family operator/(base_family lhs, const base_family &rhs)
base_family operator*(double scalar, base_family rhs)
normal_meanfield & operator/=(const normal_meanfield &rhs)
Return this approximation after elementwise division by the specified approximation's mean and log st...
Probability, optimization and sampling library.
void calc_grad(normal_meanfield &elbo_grad, M &m, Eigen::VectorXd &cont_params, int n_monte_carlo_grad, BaseRNG &rng, interface_callbacks::writer::base_writer &message_writer) const
Calculates the "blackbox" gradient with respect to both the location vector (mu) and the log-std vect...
normal_meanfield & operator=(const normal_meanfield &rhs)
Return this approximation after setting its mean vector and Cholesky factor for covariance to the val...
normal_meanfield & operator+=(double scalar)
Return this approximation after adding the specified scalar to each entry in the mean and log standar...
const Eigen::VectorXd & omega() const
Return the log standard deviation vector.
void gradient(const M &model, const Eigen::Matrix< double, Eigen::Dynamic, 1 > &x, double &f, Eigen::Matrix< double, Eigen::Dynamic, 1 > &grad_f, std::ostream *msgs=0)
Definition: util.hpp:422
normal_meanfield square() const
Return a new mean field approximation resulting from squaring the entries in the mean and log standar...
normal_meanfield & operator+=(const normal_meanfield &rhs)
Add the mean and Cholesky factor of the covariance matrix of the specified approximation to this appr...
int dimension() const
Return the dimensionality of the approximation.
base_family operator+(base_family lhs, const base_family &rhs)
normal_meanfield(const Eigen::VectorXd &mu, const Eigen::VectorXd &omega)
Construct a variational distribution with the specified mean and log standard deviation vectors...
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
const Eigen::VectorXd & mu() const
Return the mean vector.
normal_meanfield sqrt() const
Return a new mean field approximation resulting from taking the square root of the entries in the mea...
double entropy() const
Return the entropy of the approximation.
normal_meanfield(const Eigen::VectorXd &cont_params)
Construct a variational distribution with the specified mean vector and zero log standard deviation (...
void set_to_zero()
Sets the mean and log standard deviation vector for this approximation to zero.
normal_meanfield(size_t dimension)
Construct a variational distribution of the specified dimensionality with a zero mean and zero log st...
void sample(BaseRNG &rng, Eigen::VectorXd &eta) const
Assign a draw from this mean field approximation to the specified vector using the specified random n...
Variational family approximation with mean-field (diagonal covariance) multivariate normal distributi...
Eigen::VectorXd transform(const Eigen::VectorXd &eta) const
Return the transform of the sepcified vector using the Cholesky factor and mean vector.
void set_mu(const Eigen::VectorXd &mu)
Set the mean vector to the specified value.
normal_meanfield & operator*=(double scalar)
Return this approximation after multiplying by the specified scalar to each entry in the mean and log...

     [ Stan Home Page ] © 2011–2016, Stan Development Team.