1 #ifndef STAN_VARIATIONAL_NORMAL_MEANFIELD_HPP
2 #define STAN_VARIATIONAL_NORMAL_MEANFIELD_HPP
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>
27 namespace variational {
43 Eigen::VectorXd omega_;
59 : mu_(Eigen::VectorXd::Zero(dimension)),
60 omega_(Eigen::VectorXd::Zero(dimension)),
61 dimension_(dimension) {
73 omega_(Eigen::VectorXd::Zero(cont_params.size())),
74 dimension_(cont_params.size()) {
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_);
107 const Eigen::VectorXd&
mu()
const {
return mu_; }
112 const Eigen::VectorXd&
omega()
const {
return omega_; }
123 static const char*
function =
124 "stan::variational::normal_meanfield::set_mu";
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);
143 static const char*
function =
144 "stan::variational::normal_meanfield::set_omega";
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);
158 mu_ = Eigen::VectorXd::Zero(dimension_);
159 omega_ = Eigen::VectorXd::Zero(dimension_);
170 Eigen::VectorXd(omega_.array().square()));
185 Eigen::VectorXd(omega_.array().sqrt()));
200 static const char*
function =
201 "stan::variational::normal_meanfield::operator=";
202 stan::math::check_size_match(
function,
203 "Dimension of lhs", dimension_,
206 omega_ = rhs.
omega();
222 static const char*
function =
223 "stan::variational::normal_meanfield::operator+=";
224 stan::math::check_size_match(
function,
225 "Dimension of lhs", dimension_,
228 omega_ += rhs.
omega();
245 static const char*
function =
246 "stan::variational::normal_meanfield::operator/=";
247 stan::math::check_size_match(
function,
248 "Dimension of lhs", dimension_,
250 mu_.array() /= rhs.
mu().array();
251 omega_.array() /= rhs.
omega().array();
267 mu_.array() += scalar;
268 omega_.array() += scalar;
297 const Eigen::VectorXd&
mean()
const {
312 return 0.5 *
static_cast<double>(dimension_) *
313 (1.0 + stan::math::LOG_TWO_PI) + omega_.sum();
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);
336 return eta.array().cwiseProduct(omega_.array().exp()) + mu_.array();
349 template <
class BaseRNG>
350 void sample(BaseRNG& rng, Eigen::VectorXd& eta)
const {
352 for (
int d = 0; d < dimension_; ++d)
353 eta(d) = stan::math::normal_rng(0, 1, rng);
375 template <
class M,
class BaseRNG>
378 Eigen::VectorXd& cont_params,
379 int n_monte_carlo_grad,
383 static const char*
function =
384 "stan::variational::normal_meanfield::calc_grad";
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());
393 Eigen::VectorXd mu_grad = Eigen::VectorXd::Zero(dimension_);
394 Eigen::VectorXd omega_grad = Eigen::VectorXd::Zero(dimension_);
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_);
401 static const int n_retries = 10;
402 for (
int i = 0, n_monte_carlo_drop = 0; i < n_monte_carlo_grad; ) {
404 for (
int d = 0; d < dimension_; ++d)
405 eta(d) = stan::math::normal_rng(0, 1, rng);
408 std::stringstream 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());
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);
428 mu_grad /=
static_cast<double>(n_monte_carlo_grad);
429 omega_grad /=
static_cast<double>(n_monte_carlo_grad);
432 = omega_grad.array().cwiseProduct(omega_.array().exp());
434 omega_grad.array() += 1.0;
436 elbo_grad.
set_mu(mu_grad);
479 return rhs += scalar;
492 return rhs *= scalar;
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)
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.
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...