Test (adaptive) Metropolis-Hastings with diagonally elongated distribution

[dawsoneliasen]

WIP. @bangerth could you help me out with the #FIXMEs here?

[bangerth]

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?

[bangerth]

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.

[bangerth]

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

[bangerth]

double log_likelihood (const SampleType &x)
{
  double mu[2] = {0, 0};
  // double cov[2][2] = {{10, 0}, {0, 1}};

  double cov_inv[2][2] = {{1./100, 0},
                          {0, 1}};
  double phi = 3.1415/6;
  double rotation_matrix[2][2] = {{cos(phi), sin(phi)},
                                  {-sin(phi), cos(phi)}};
  double rotated_cov_inv = rotation_matrix * cov_inv * rotation_matrix.transpose();
  
  
  double p =  -0.5 * ((x[0]-mu[0])*cov_inv[0][0]*(x[0]-mu[0]) +
                      (x[0]-mu[0])*cov_inv[0][1]*(x[1]-mu[1]) +
                      (x[1]-mu[1])*cov_inv[1][0]*(x[0]-mu[0]) +
                      (x[1]-mu[1])*cov_inv[1][1]*(x[1]-mu[1]));
  std::cout << x[0] << " " << x[1] << " " << p << std::endl;
  return p;
}


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

[dawsoneliasen]

@bangerth PTAL. Note that I replaced the Cholesky implementation with something much simpler (took better advantage of the fact that I'm assuming A is 2x2 matrix). I think the covariance is correct but I'm not positive -- recall that we took an elongated distribution and "rotated" it by a kind of arbitrary factor. Appreciate any insight here.

[bangerth]

Can you compute what the correct covariance matrix is from the one you put in, and compare to what you get out?

[dawsoneliasen]

@bangerth Yes, good idea. My results are close to the results of the test:

7.75    -3.90
-3.90    3.25

[dawsoneliasen]

@bangerth And, when in increase the number of samples in the test by a factor of ten, the test output gets a lot closer to these values.

[bangerth]

OK, so you think this is good to merge?