From: Wolfgang Bangerth Date: Thu, 2 Jan 2020 01:06:46 +0000 (-0700) Subject: Fix a bug in the Metropolis-Hastings solver. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b11158d4bf42d0faaf419974deabbf07d6ba5a98;p=code-gallery.git Fix a bug in the Metropolis-Hastings solver. --- diff --git a/MCMC-Laplace/mcmc-laplace.cc b/MCMC-Laplace/mcmc-laplace.cc index 813f048..47417a0 100644 --- a/MCMC-Laplace/mcmc-laplace.cc +++ b/MCMC-Laplace/mcmc-laplace.cc @@ -74,6 +74,8 @@ namespace ForwardSimulator virtual Vector evaluate(const Vector &coefficients) = 0; }; + + template class PoissonSolver : public Interface { @@ -499,12 +501,26 @@ namespace LogPrior // to multiply the existing sample entries by are close to one. And // because the exponential of a number is always positive, we never // get negative samples this way.) +// +// But the Metropolis-Hastings sampler doesn't just need a perturbed +// sample $y$ location given the current sample location $x$. It also +// needs to know the ratio of the probability of reaching $y$ from +// $x$, divided by the probability of reaching $x$ from $y$. If we +// were to use a symmetric proposal distribution (e.g., a Gaussian +// distribution centered at $x$ with a width independent of $x$), then +// these two probabilities would be the same, and the ratio one. But +// that's not the case for the Gaussian in log space. It's not +// terribly difficult to verify that in that case, for a single +// component the ratio of these probabilities is $y_i/x_i$, and +// consequently for all components of the vector together, the +// probability is the product of these ratios. namespace ProposalGenerator { class Interface { public: - virtual Vector + virtual + std::pair,double> perturb(const Vector ¤t_sample) const = 0; }; @@ -514,7 +530,9 @@ namespace ProposalGenerator public: LogGaussian(const unsigned int random_seed, const double log_sigma); - virtual Vector perturb(const Vector ¤t_sample) const; + virtual + std::pair,double> + perturb(const Vector ¤t_sample) const; private: const double log_sigma; @@ -522,6 +540,7 @@ namespace ProposalGenerator }; + LogGaussian::LogGaussian(const unsigned int random_seed, const double log_sigma) : log_sigma(log_sigma) @@ -529,15 +548,21 @@ namespace ProposalGenerator random_number_generator.seed(random_seed); } - Vector + + std::pair,double> LogGaussian::perturb(const Vector ¤t_sample) const { Vector new_sample = current_sample; + double product_of_ratios = 1; for (auto &x : new_sample) - x *= std::exp( - std::normal_distribution<>(0, log_sigma)(random_number_generator)); + { + const double rnd = std::normal_distribution<>(0, log_sigma)(random_number_generator); + const double exp_rnd = std::exp(rnd); + x *= exp_rnd; + product_of_ratios *= exp_rnd; + } - return new_sample; + return {new_sample, product_of_ratios}; } } // namespace ProposalGenerator @@ -624,14 +649,22 @@ namespace Sampler for (unsigned int k = 1; k < n_samples; ++k, ++sample_number) { - const Vector trial_sample = - proposal_generator.perturb(current_sample); + std::pair,double> + perturbation = proposal_generator.perturb(current_sample); + const Vector trial_sample = std::move (perturbation.first); + const double perturbation_probability_ratio = perturbation.second; + const double trial_log_posterior = (likelihood.log_likelihood(simulator.evaluate(trial_sample)) + prior.log_prior(trial_sample)); - if ((trial_log_posterior > current_log_posterior) || - (std::exp(trial_log_posterior - current_log_posterior) >= + if ((trial_log_posterior + std::log(perturbation_probability_ratio) + >= + current_log_posterior) + || + (std::exp(trial_log_posterior - current_log_posterior) + * perturbation_probability_ratio + >= uniform_distribution(random_number_generator))) { current_sample = trial_sample;