]> https://gitweb.dealii.org/ - code-gallery.git/commitdiff
Fix a bug in the Metropolis-Hastings solver.
authorWolfgang Bangerth <bangerth@colostate.edu>
Thu, 2 Jan 2020 01:06:46 +0000 (18:06 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 16 Jan 2020 22:02:43 +0000 (15:02 -0700)
MCMC-Laplace/mcmc-laplace.cc

index 813f048cd39ac14da1d09999f270ca5850a874a3..47417a0e30f0fabc1516d618c42cf19194144afc 100644 (file)
@@ -74,6 +74,8 @@ namespace ForwardSimulator
     virtual Vector<double> evaluate(const Vector<double> &coefficients) = 0;
   };
 
+
+
   template <int dim>
   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<double>
+    virtual
+    std::pair<Vector<double>,double>
     perturb(const Vector<double> &current_sample) const = 0;
   };
 
@@ -514,7 +530,9 @@ namespace ProposalGenerator
   public:
     LogGaussian(const unsigned int random_seed, const double log_sigma);
 
-    virtual Vector<double> perturb(const Vector<double> &current_sample) const;
+    virtual
+    std::pair<Vector<double>,double>
+    perturb(const Vector<double> &current_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<double>
+
+  std::pair<Vector<double>,double>
   LogGaussian::perturb(const Vector<double> &current_sample) const
   {
     Vector<double> 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<double> trial_sample =
-          proposal_generator.perturb(current_sample);
+        std::pair<Vector<double>,double>
+          perturbation = proposal_generator.perturb(current_sample);
+        const Vector<double> 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;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.