Skip to content

Test (adaptive) Metropolis-Hastings with diagonally elongated distribution #148

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dawsoneliasen wants to merge 17 commits into
base: master
Choose a base branch
from
Open

Test (adaptive) Metropolis-Hastings with diagonally elongated distribution #148

Changes from 1 commit
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
624df6e
create new met-has test
dawsoneliasen Dec 7, 2020
1743db4
implement likelihood
dawsoneliasen Dec 15, 2020
2fb59e2
remove comment
dawsoneliasen Dec 16, 2020
2201aca
add adaptive test
dawsoneliasen Dec 16, 2020
1f74e61
use inverse of covariance matrix
dawsoneliasen Dec 19, 2020
79e9482
add output
dawsoneliasen Jan 21, 2021
bb0220f
fix perturb, rotate cov
dawsoneliasen Jan 28, 2021
34b6ad6
copy over changes from 12
dawsoneliasen Jan 28, 2021
c0f9de9
fix compiler errors
dawsoneliasen Jan 28, 2021
3547070
use ublas matrices
dawsoneliasen Feb 4, 2021
8abc056
use diagonal matrix, not vector
dawsoneliasen Feb 4, 2021
5040a9e
remove comments
dawsoneliasen Feb 4, 2021
463a8e1
use static RNG and scale the normal vector
dawsoneliasen Feb 9, 2021
c175fa2
add covariance to output
dawsoneliasen Feb 9, 2021
ad199e2
fix cholesky
dawsoneliasen Feb 9, 2021
9220100
add output files
dawsoneliasen Feb 9, 2021
69fba75
update test 12 output
dawsoneliasen Feb 9, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 63 additions & 0 deletions tests/metropolis_hasting_producer_12.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2020 by the SampleFlow authors.
//
// This file is part of the SampleFlow library.
//
// The SampleFlow library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of SampleFlow.
//
// ---------------------------------------------------------------------


// Test the Metropolis-Hastings producer with a diagonally elongated
// distribution.

#include <iostream>
#include <sampleflow/producers/metropolis_hastings.h>
#include <sampleflow/filters/conversion.h>
#include <sampleflow/consumers/mean_value.h>
#include <valarray>
#include <random>
#include <cmath>


using SampleType = std::valarray<double>;


double log_likelihood (const SampleType &x)
{
std::valarray<double> mu = {0, 0};
std::valarray<double> cov = {{1, 1.5}, {1.5, 3}};

@bangerth bangerth Dec 10, 2020

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This won't compile. But you can just do

  double mu[2] = {0,0};
  double cov[2][2] = {{1,1.5}, {1.5,3}};

Have you checked whether this really corresponds to positive definite matrix?

// FIXME: what's the correct way to do these matrix operations?
return -0.5 * (ln(|cov|) + (x - mu)^T @ cov^-1 @ (x - mu) + k * ln(2 * pi));

@bangerth bangerth Dec 10, 2020

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

First, the log_likelihood doesn't have to be normalized, so you can remove the ln(|cov|) factor along with the ln(2*pi) one.

Second, just multiply it out:

Suggested change
return -0.5 * (ln(|cov|) + (x - mu)^T @ cov^-1 @ (x - mu) + k * ln(2 * pi));
return -0.5 * ( (x[0]-mu[0])*cov[0][0]*(x[0]-mu[0]) +
(x[0]-mu[0])*cov[0][1]*(x[1]-mu[1]) +
(x[1]-mu[1])*cov[1][0]*(x[0]-mu[0]) +
(x[1]-mu[1])*cov[1][1]*(x[1]-mu[1]) +

That's not elegant, but it'll do for this test.

}


std::pair<SampleType,double> perturb (const SampleType &x)
{
// FIXME: how to sample MVN(x, cov)?

@bangerth bangerth Dec 10, 2020

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

as discussed today, you don't have to :-)

static std::mt19937 rng;
std::normal_distribution<double> distribution(0, 1);
SampleType y = x;
for (auto &el : y)
y += distribution(rng);
return {y, 1.0};
}


int main ()
{
SampleFlow::Producers::MetropolisHastings<SampleType> mh_sampler;
SampleFlow::Producers::CovarianceMatrix<SampleType> cov_matrix;
SampleFlow::Consumers::MeanValue<SampleType> mean_value;
cov_matrix.connect_to_producer(mh_sampler);
mean_value.connect_to_producer(mh_sampler);
// Sample, starting at an asymmetric point, and creating 100,000 samples
mh_sampler.sample({5, -1}, &log_likelihood, &perturb, 100000);
std::cout << "Mean value = " << mean_value.get() << std::endl;
}