virtual Vector<double> evaluate(const Vector<double> &coefficients) = 0;
};
+
+
template <int dim>
class PoissonSolver : public Interface
{
// 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> ¤t_sample) const = 0;
};
public:
LogGaussian(const unsigned int random_seed, const double log_sigma);
- virtual Vector<double> perturb(const Vector<double> ¤t_sample) const;
+ virtual
+ std::pair<Vector<double>,double>
+ perturb(const Vector<double> ¤t_sample) const;
private:
const double log_sigma;
};
+
LogGaussian::LogGaussian(const unsigned int random_seed,
const double log_sigma)
: log_sigma(log_sigma)
random_number_generator.seed(random_seed);
}
- Vector<double>
+
+ std::pair<Vector<double>,double>
LogGaussian::perturb(const Vector<double> ¤t_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
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;