Stan  2.14.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.hpp>
8 #include <algorithm>
9 #include <ostream>
10 #include <vector>
11 
12 namespace stan {
13 
14  namespace variational {
15 
20  class normal_meanfield : public base_family {
21  private:
25  Eigen::VectorXd mu_;
26 
30  Eigen::VectorXd omega_;
31 
35  const int dimension_;
36 
37  public:
45  explicit normal_meanfield(size_t dimension)
46  : mu_(Eigen::VectorXd::Zero(dimension)),
47  omega_(Eigen::VectorXd::Zero(dimension)),
48  dimension_(dimension) {
49  }
50 
58  explicit normal_meanfield(const Eigen::VectorXd& cont_params)
59  : mu_(cont_params),
60  omega_(Eigen::VectorXd::Zero(cont_params.size())),
61  dimension_(cont_params.size()) {
62  }
63 
74  normal_meanfield(const Eigen::VectorXd& mu,
75  const Eigen::VectorXd& omega)
76  : mu_(mu), omega_(omega), dimension_(mu.size()) {
77  static const char* function =
78  "stan::variational::normal_meanfield";
79  stan::math::check_size_match(function,
80  "Dimension of mean vector", dimension_,
81  "Dimension of log std vector", omega_.size() );
82  stan::math::check_not_nan(function, "Mean vector", mu_);
83  stan::math::check_not_nan(function, "Log std vector", omega_);
84  }
85 
89  int dimension() const { return dimension_; }
90 
94  const Eigen::VectorXd& mu() const { return mu_; }
95 
99  const Eigen::VectorXd& omega() const { return omega_; }
100 
109  void set_mu(const Eigen::VectorXd& mu) {
110  static const char* function =
111  "stan::variational::normal_meanfield::set_mu";
112 
113  stan::math::check_size_match(function,
114  "Dimension of input vector", mu.size(),
115  "Dimension of current vector", dimension_);
116  stan::math::check_not_nan(function, "Input vector", mu);
117  mu_ = mu;
118  }
119 
129  void set_omega(const Eigen::VectorXd& omega) {
130  static const char* function =
131  "stan::variational::normal_meanfield::set_omega";
132 
133  stan::math::check_size_match(function,
134  "Dimension of input vector", omega.size(),
135  "Dimension of current vector", dimension_);
136  stan::math::check_not_nan(function, "Input vector", omega);
137  omega_ = omega;
138  }
139 
144  void set_to_zero() {
145  mu_ = Eigen::VectorXd::Zero(dimension_);
146  omega_ = Eigen::VectorXd::Zero(dimension_);
147  }
148 
156  return normal_meanfield(Eigen::VectorXd(mu_.array().square()),
157  Eigen::VectorXd(omega_.array().square()));
158  }
159 
171  return normal_meanfield(Eigen::VectorXd(mu_.array().sqrt()),
172  Eigen::VectorXd(omega_.array().sqrt()));
173  }
174 
187  static const char* function =
188  "stan::variational::normal_meanfield::operator=";
189  stan::math::check_size_match(function,
190  "Dimension of lhs", dimension_,
191  "Dimension of rhs", rhs.dimension());
192  mu_ = rhs.mu();
193  omega_ = rhs.omega();
194  return *this;
195  }
196 
209  static const char* function =
210  "stan::variational::normal_meanfield::operator+=";
211  stan::math::check_size_match(function,
212  "Dimension of lhs", dimension_,
213  "Dimension of rhs", rhs.dimension());
214  mu_ += rhs.mu();
215  omega_ += rhs.omega();
216  return *this;
217  }
218 
232  static const char* function =
233  "stan::variational::normal_meanfield::operator/=";
234  stan::math::check_size_match(function,
235  "Dimension of lhs", dimension_,
236  "Dimension of rhs", rhs.dimension());
237  mu_.array() /= rhs.mu().array();
238  omega_.array() /= rhs.omega().array();
239  return *this;
240  }
241 
253  normal_meanfield& operator+=(double scalar) {
254  mu_.array() += scalar;
255  omega_.array() += scalar;
256  return *this;
257  }
258 
271  normal_meanfield& operator*=(double scalar) {
272  mu_ *= scalar;
273  omega_ *= scalar;
274  return *this;
275  }
276 
284  const Eigen::VectorXd& mean() const {
285  return mu();
286  }
287 
298  double entropy() const {
299  return 0.5 * static_cast<double>(dimension_) *
300  (1.0 + stan::math::LOG_TWO_PI) + omega_.sum();
301  }
302 
315  Eigen::VectorXd transform(const Eigen::VectorXd& eta) const {
316  static const char* function =
317  "stan::variational::normal_meanfield::transform";
318  stan::math::check_size_match(function,
319  "Dimension of mean vector", dimension_,
320  "Dimension of input vector", eta.size() );
321  stan::math::check_not_nan(function, "Input vector", eta);
322  // exp(omega) * eta + mu
323  return eta.array().cwiseProduct(omega_.array().exp()) + mu_.array();
324  }
325 
336  template <class BaseRNG>
337  void sample(BaseRNG& rng, Eigen::VectorXd& eta) const {
338  // Draw from standard normal and transform to real-coordinate space
339  for (int d = 0; d < dimension_; ++d)
340  eta(d) = stan::math::normal_rng(0, 1, rng);
341  eta = transform(eta);
342  }
343 
362  template <class M, class BaseRNG>
363  void calc_grad(normal_meanfield& elbo_grad,
364  M& m,
365  Eigen::VectorXd& cont_params,
366  int n_monte_carlo_grad,
367  BaseRNG& rng,
369  const {
370  static const char* function =
371  "stan::variational::normal_meanfield::calc_grad";
372 
373  stan::math::check_size_match(function,
374  "Dimension of elbo_grad", elbo_grad.dimension(),
375  "Dimension of variational q", dimension_);
376  stan::math::check_size_match(function,
377  "Dimension of variational q", dimension_,
378  "Dimension of variables in model", cont_params.size());
379 
380  Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
381  Eigen::VectorXd omega_grad = Eigen::VectorXd::Zero(dimension_);
382  double tmp_lp = 0.0;
383  Eigen::VectorXd tmp_mu_grad = Eigen::VectorXd::Zero(dimension_);
384  Eigen::VectorXd eta = Eigen::VectorXd::Zero(dimension_);
385  Eigen::VectorXd zeta = Eigen::VectorXd::Zero(dimension_);
386 
387  // Naive Monte Carlo integration
388  static const int n_retries = 10;
389  for (int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
390  // Draw from standard normal and transform to real-coordinate space
391  for (int d = 0; d < dimension_; ++d)
392  eta(d) = stan::math::normal_rng(0, 1, rng);
393  zeta = transform(eta);
394  try {
395  std::stringstream ss;
396  stan::model::gradient(m, zeta, tmp_lp, tmp_mu_grad, &ss);
397  if (ss.str().length() > 0)
398  message_writer(ss.str());
399  stan::math::check_finite(function, "Gradient of mu", tmp_mu_grad);
400  mu_grad += tmp_mu_grad;
401  omega_grad.array() += tmp_mu_grad.array().cwiseProduct(eta.array());
402  ++i;
403  } catch (const std::exception& e) {
404  ++n_monte_carlo_drop;
405  if (n_monte_carlo_drop >= n_retries * n_monte_carlo_grad) {
406  const char* name = "The number of dropped evaluations";
407  const char* msg1 = "has reached its maximum amount (";
408  int y = n_retries * n_monte_carlo_grad;
409  const char* msg2 = "). Your model may be either severely "
410  "ill-conditioned or misspecified.";
411  stan::math::domain_error(function, name, y, msg1, msg2);
412  }
413  }
414  }
415  mu_grad /= static_cast<double>(n_monte_carlo_grad);
416  omega_grad /= static_cast<double>(n_monte_carlo_grad);
417 
418  omega_grad.array()
419  = omega_grad.array().cwiseProduct(omega_.array().exp());
420 
421  omega_grad.array() += 1.0; // add entropy gradient (unit)
422 
423  elbo_grad.set_mu(mu_grad);
424  elbo_grad.set_omega(omega_grad);
425  }
426  };
427 
438  const normal_meanfield& rhs) {
439  return lhs += rhs;
440  }
441 
452  const normal_meanfield& rhs) {
453  return lhs /= rhs;
454  }
455 
466  return rhs += scalar;
467  }
468 
479  return rhs *= scalar;
480  }
481 
482  }
483 }
484 #endif
void set_omega(const Eigen::VectorXd &omega)
Set the log standard deviation vector to the specified value.
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&#39;s mean and log st...
Probability, optimization and sampling library.
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...
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...
normal_meanfield sqrt() const
Return a new mean field approximation resulting from taking the square root of the entries in the mea...
const Eigen::VectorXd & omega() const
Return the log standard deviation vector.
double entropy() const
Return the entropy of the approximation.
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: gradient.hpp:14
normal_meanfield & operator+=(const normal_meanfield &rhs)
Add the mean and Cholesky factor of the covariance matrix of the specified approximation to this appr...
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...
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...
normal_meanfield square() const
Return a new mean field approximation resulting from squaring the entries in the mean and log standar...
base_writer is an abstract base class defining the interface for Stan writer callbacks.
Definition: base_writer.hpp:20
normal_meanfield(const Eigen::VectorXd &cont_params)
Construct a variational distribution with the specified mean vector and zero log standard deviation (...
Eigen::VectorXd transform(const Eigen::VectorXd &eta) const
Return the transform of the sepcified vector using the Cholesky factor and mean vector.
const Eigen::VectorXd & mean() const
Returns the mean vector for this approximation.
void set_to_zero()
Sets the mean and log standard deviation vector for this approximation to zero.
int dimension() const
Return the dimensionality of the approximation.
normal_meanfield(size_t dimension)
Construct a variational distribution of the specified dimensionality with a zero mean and zero log st...
const Eigen::VectorXd & mu() const
Return the mean vector.
Variational family approximation with mean-field (diagonal covariance) multivariate normal distributi...
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.