From 624df6ed73a9d39f02f1fd594fee3388b8537146 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Mon, 7 Dec 2020 09:37:41 -0700 Subject: [PATCH 01/17] create new met-has test --- tests/metropolis_hasting_producer_12.cc | 63 +++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 tests/metropolis_hasting_producer_12.cc diff --git a/tests/metropolis_hasting_producer_12.cc b/tests/metropolis_hasting_producer_12.cc new file mode 100644 index 00000000..fd653a87 --- /dev/null +++ b/tests/metropolis_hasting_producer_12.cc @@ -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 +#include +#include +#include +#include +#include +#include + + +using SampleType = std::valarray; + + +double log_likelihood (const SampleType &x) +{ + std::valarray mu = {0, 0}; + std::valarray cov = {{1, 1.5}, {1.5, 3}}; + // 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)); +} + + +std::pair perturb (const SampleType &x) +{ + // FIXME: how to sample MVN(x, cov)? + static std::mt19937 rng; + std::normal_distribution distribution(0, 1); + SampleType y = x; + for (auto &el : y) + y += distribution(rng); + return {y, 1.0}; +} + + +int main () +{ + SampleFlow::Producers::MetropolisHastings mh_sampler; + SampleFlow::Producers::CovarianceMatrix cov_matrix; + SampleFlow::Consumers::MeanValue 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; +} From 1743db49dc689ee3b5d67cb125f167e6e60beb8f Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Tue, 15 Dec 2020 08:46:48 -0700 Subject: [PATCH 02/17] implement likelihood --- tests/metropolis_hasting_producer_12.cc | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/tests/metropolis_hasting_producer_12.cc b/tests/metropolis_hasting_producer_12.cc index fd653a87..74f3ed8c 100644 --- a/tests/metropolis_hasting_producer_12.cc +++ b/tests/metropolis_hasting_producer_12.cc @@ -15,11 +15,13 @@ // Test the Metropolis-Hastings producer with a diagonally elongated -// distribution. +// distribution; here, we use a MVN with mean [0, 0] and covariance +// [[10, 0], [0, 1]]. #include #include #include +#include #include #include #include @@ -31,16 +33,18 @@ using SampleType = std::valarray; double log_likelihood (const SampleType &x) { - std::valarray mu = {0, 0}; - std::valarray cov = {{1, 1.5}, {1.5, 3}}; - // 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)); + double mu[2] = {0, 0}; + double cov[2][2] = {{10, 0}, {0, 1}}; + // 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])); } std::pair perturb (const SampleType &x) { - // FIXME: how to sample MVN(x, cov)? static std::mt19937 rng; std::normal_distribution distribution(0, 1); SampleType y = x; @@ -53,11 +57,12 @@ std::pair perturb (const SampleType &x) int main () { SampleFlow::Producers::MetropolisHastings mh_sampler; - SampleFlow::Producers::CovarianceMatrix cov_matrix; + SampleFlow::Consumers::CovarianceMatrix cov_matrix; SampleFlow::Consumers::MeanValue 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; + std::cout << "Mean value = " << mean_value.get()[0] << std::endl; + std::cout << "Mean value = " << mean_value.get()[1] << std::endl; } From 2fb59e25313f609b3d007722da8015a3d98e3d6f Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Wed, 16 Dec 2020 09:52:14 -0700 Subject: [PATCH 03/17] remove comment --- tests/metropolis_hasting_producer_12.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/metropolis_hasting_producer_12.cc b/tests/metropolis_hasting_producer_12.cc index 74f3ed8c..16ca0df7 100644 --- a/tests/metropolis_hasting_producer_12.cc +++ b/tests/metropolis_hasting_producer_12.cc @@ -35,7 +35,6 @@ double log_likelihood (const SampleType &x) { double mu[2] = {0, 0}; double cov[2][2] = {{10, 0}, {0, 1}}; - // 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]) + From 2201aca6acbe588eebee0b2c212da7ee2cca9277 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Wed, 16 Dec 2020 15:05:48 -0700 Subject: [PATCH 04/17] add adaptive test --- tests/metropolis_hasting_producer_13.cc | 107 ++++++++++++++++++++++++ 1 file changed, 107 insertions(+) create mode 100644 tests/metropolis_hasting_producer_13.cc diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc new file mode 100644 index 00000000..d0d33356 --- /dev/null +++ b/tests/metropolis_hasting_producer_13.cc @@ -0,0 +1,107 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// Copy of metropolis_hasting_producer_12, but with adaptive sampling. +// Test the Metropolis-Hastings producer with a diagonally elongated +// distribution; here, we use a MVN with mean [0, 0] and covariance +// [[10, 0], [0, 1]]. + +#include +#include +#include +#include +#include +#include +#include +#include + + +using SampleType = std::valarray; + + +double log_likelihood (const SampleType &x) +{ + double mu[2] = {0, 0}; + double cov[2][2] = {{10, 0}, {0, 1}}; + // 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])); +} + + +std::pair cholesky(A) +{ + double L[A.size()][A[0].size()]; + double D[A[0].size()]; + for(i = 0; i < A.size(); ++i) + { + for (j = 0; j < A[0].size(); ++j) + { + double partial_sum = 0; + for (k = 0; k < (j - 1); ++k) + partial_sum += L[j][k] ** 2 * D[k]; + D[j] = A[j][j] - partial_sum; + if (i > j) + { + for (k = 0; k < (j - 1); ++k) + partial_sum += L[i][k] * L[j][k] * D[k]; + L[i][j] = 1 / D[j] * (A[i][j] - partial_sum); + } + else if (i == j) + L[i][j] = 1; + else + L[i][j] = 0; + } + } + return {L, D}; +} + + +std::pair perturb (const SampleType &x, const double[] &cov) +{ + std::pair decomp = cholesky(cov); + double L[][] = decomp.first; + double D[] = decomp.second; + static std::mt19937 rng; + SampleType y = x; + for (i = 0; i < y.size(); ++i) + { + std::normal_distribution distribution(x[i], D[i]); + y[i] = L[i][j] * distribution(rng) * L[j][i]; + } + return {y, 1.0}; +} + + +int main () +{ + SampleFlow::Producers::MetropolisHastings mh_sampler; + SampleFlow::Consumers::CovarianceMatrix cov_matrix; + SampleFlow::Consumers::MeanValue 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, + [](const SampleType &x){return perturb(x, cov_matrix.get())}, + 100000 + ); + std::cout << "Mean value = " << mean_value.get()[0] << std::endl; + std::cout << "Mean value = " << mean_value.get()[1] << std::endl; +} From 1f74e61ba7a251542ca20002a7776a5780f8ed67 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Sat, 19 Dec 2020 09:26:47 -0700 Subject: [PATCH 05/17] use inverse of covariance matrix --- tests/metropolis_hasting_producer_12.cc | 11 ++++++----- tests/metropolis_hasting_producer_13.cc | 11 +++++++---- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/tests/metropolis_hasting_producer_12.cc b/tests/metropolis_hasting_producer_12.cc index 16ca0df7..1549f284 100644 --- a/tests/metropolis_hasting_producer_12.cc +++ b/tests/metropolis_hasting_producer_12.cc @@ -34,11 +34,12 @@ using SampleType = std::valarray; double log_likelihood (const SampleType &x) { double mu[2] = {0, 0}; - double cov[2][2] = {{10, 0}, {0, 1}}; - 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])); + // double cov[2][2] = {{10, 0}, {0, 1}}; + double cov_inv[2][2] = {{1/10, 0},{0, 1}}; + return -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])); } diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index d0d33356..954d2d8d 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -75,14 +75,17 @@ std::pair cholesky(A) std::pair perturb (const SampleType &x, const double[] &cov) { std::pair decomp = cholesky(cov); - double L[][] = decomp.first; - double D[] = decomp.second; + double L[x.size()][x.size()] = decomp.first; + double D[x.size()] = decomp.second; static std::mt19937 rng; SampleType y = x; for (i = 0; i < y.size(); ++i) { - std::normal_distribution distribution(x[i], D[i]); - y[i] = L[i][j] * distribution(rng) * L[j][i]; + std::normal_distribution distribution(x[i], 1 ) * D[i]; + double perturbation[2] = distribution(rng); + // xt = xk + l * sqrt(D) * n // n is a vector whose elements are normally distrubitued with variance 1 + // + y[i] = L[i][j] * distribution(rng); } return {y, 1.0}; } From 79e94825a956428a5deb7b312fa9e252da1ae846 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 21 Jan 2021 16:31:02 -0700 Subject: [PATCH 06/17] add output --- tests/metropolis_hasting_producer_12.cc | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/tests/metropolis_hasting_producer_12.cc b/tests/metropolis_hasting_producer_12.cc index 1549f284..79a70f8b 100644 --- a/tests/metropolis_hasting_producer_12.cc +++ b/tests/metropolis_hasting_producer_12.cc @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -36,10 +37,16 @@ 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/10, 0},{0, 1}}; - return -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])); + double dev[2] = {x[0]- mu[0], x[1] - mu[1]}; + double dev_covinv[2] = {dev[0] * cov_inv[0][0] + dev[1] * cov_inv[0][1], + dev[0] * cov_inv[1][0] + dev[1] * cov_inv[1][1]}; + double dev_covinv_dev = dev_covinv[0] * dev[0] + dev_covinv[1] * dev[1]; + // return -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] << " " << -0.5 * dev_covinv_dev << std::endl; + return -0.5 * dev_covinv_dev; } @@ -61,6 +68,8 @@ int main () SampleFlow::Consumers::MeanValue mean_value; cov_matrix.connect_to_producer(mh_sampler); mean_value.connect_to_producer(mh_sampler); + // SampleFlow::Consumers::StreamOutput output(std::cout); + // output.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()[0] << std::endl; From bb0220f242d32338f3869193f2e6dcf1f2c10389 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 28 Jan 2021 14:34:06 -0700 Subject: [PATCH 07/17] fix perturb, rotate cov --- tests/metropolis_hasting_producer_12.cc | 48 +++++++++++++++++-------- 1 file changed, 34 insertions(+), 14 deletions(-) diff --git a/tests/metropolis_hasting_producer_12.cc b/tests/metropolis_hasting_producer_12.cc index 79a70f8b..6b7511aa 100644 --- a/tests/metropolis_hasting_producer_12.cc +++ b/tests/metropolis_hasting_producer_12.cc @@ -15,8 +15,7 @@ // Test the Metropolis-Hastings producer with a diagonally elongated -// distribution; here, we use a MVN with mean [0, 0] and covariance -// [[10, 0], [0, 1]]. +// distribution. #include #include @@ -36,16 +35,39 @@ 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/10, 0},{0, 1}}; - double dev[2] = {x[0]- mu[0], x[1] - mu[1]}; - double dev_covinv[2] = {dev[0] * cov_inv[0][0] + dev[1] * cov_inv[0][1], - dev[0] * cov_inv[1][0] + dev[1] * cov_inv[1][1]}; + double cov_inv[2][2] = {{1.0/10.0, 0},{0, 1}}; + + double phi = 3.1415 / 6; + double rotation_matrix[2][2] = {{cos(phi), sin(phi)}, + {-sin(phi), cos(phi)} + }; + double rotation_matrix_transpose[2][2] = {{cos(phi), -sin(phi)}, + {sin(phi), cos(phi)} + }; + double RC[2][2] = {{ + rotation_matrix[0][0] *cov_inv[0][0] + rotation_matrix[0][1] *cov_inv[1][0], + rotation_matrix[0][0] *cov_inv[0][1] + rotation_matrix[0][1] *cov_inv[1][1] + }, + { + rotation_matrix[1][0] *cov_inv[0][0] + rotation_matrix[1][1] *cov_inv[1][0], + rotation_matrix[1][0] *cov_inv[0][1] + rotation_matrix[1][1] *cov_inv[1][1] + } + }; + double rotated_cov_inv[2][2] = {{ + RC[0][0] *rotation_matrix_transpose[0][0] + RC[0][1] *rotation_matrix_transpose[1][0], + RC[0][0] *rotation_matrix_transpose[0][1] + RC[0][1] *rotation_matrix_transpose[1][1] + }, + { + RC[1][0] *rotation_matrix_transpose[0][0] + RC[1][1] *rotation_matrix_transpose[1][0], + RC[1][0] *rotation_matrix_transpose[0][1] + RC[1][1] *rotation_matrix_transpose[1][1] + } + }; + + double dev[2] = {x[0] - mu[0], x[1] - mu[1]}; + double dev_covinv[2] = {dev[0] *rotated_cov_inv[0][0] + dev[1] *rotated_cov_inv[0][1], + dev[0] *rotated_cov_inv[1][0] + dev[1] *rotated_cov_inv[1][1] + }; double dev_covinv_dev = dev_covinv[0] * dev[0] + dev_covinv[1] * dev[1]; - // return -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] << " " << -0.5 * dev_covinv_dev << std::endl; return -0.5 * dev_covinv_dev; } @@ -56,7 +78,7 @@ std::pair perturb (const SampleType &x) std::normal_distribution distribution(0, 1); SampleType y = x; for (auto &el : y) - y += distribution(rng); + el += distribution(rng); return {y, 1.0}; } @@ -68,8 +90,6 @@ int main () SampleFlow::Consumers::MeanValue mean_value; cov_matrix.connect_to_producer(mh_sampler); mean_value.connect_to_producer(mh_sampler); - // SampleFlow::Consumers::StreamOutput output(std::cout); - // output.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()[0] << std::endl; From 34b6ad69e18c5fcc42b1a78b64207900fd025170 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 28 Jan 2021 14:34:15 -0700 Subject: [PATCH 08/17] copy over changes from 12 --- tests/metropolis_hasting_producer_13.cc | 99 ++++++++++++++++--------- 1 file changed, 66 insertions(+), 33 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index 954d2d8d..9828fcc5 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -16,8 +16,7 @@ // Copy of metropolis_hasting_producer_12, but with adaptive sampling. // Test the Metropolis-Hastings producer with a diagonally elongated -// distribution; here, we use a MVN with mean [0, 0] and covariance -// [[10, 0], [0, 1]]. +// distribution. #include #include @@ -35,30 +34,59 @@ using SampleType = std::valarray; double log_likelihood (const SampleType &x) { double mu[2] = {0, 0}; - double cov[2][2] = {{10, 0}, {0, 1}}; - // 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])); + // double cov[2][2] = {{10, 0}, {0, 1}}; + double cov_inv[2][2] = {{1.0/10.0, 0},{0, 1}}; + + double phi = 3.1415 / 6; + double rotation_matrix[2][2] = {{cos(phi), sin(phi)}, + {-sin(phi), cos(phi)} + }; + double rotation_matrix_transpose[2][2] = {{cos(phi), -sin(phi)}, + {sin(phi), cos(phi)} + }; + double RC[2][2] = {{ + rotation_matrix[0][0] *cov_inv[0][0] + rotation_matrix[0][1] *cov_inv[1][0], + rotation_matrix[0][0] *cov_inv[0][1] + rotation_matrix[0][1] *cov_inv[1][1] + }, + { + rotation_matrix[1][0] *cov_inv[0][0] + rotation_matrix[1][1] *cov_inv[1][0], + rotation_matrix[1][0] *cov_inv[0][1] + rotation_matrix[1][1] *cov_inv[1][1] + } + }; + double rotated_cov_inv[2][2] = {{ + RC[0][0] *rotation_matrix_transpose[0][0] + RC[0][1] *rotation_matrix_transpose[1][0], + RC[0][0] *rotation_matrix_transpose[0][1] + RC[0][1] *rotation_matrix_transpose[1][1] + }, + { + RC[1][0] *rotation_matrix_transpose[0][0] + RC[1][1] *rotation_matrix_transpose[1][0], + RC[1][0] *rotation_matrix_transpose[0][1] + RC[1][1] *rotation_matrix_transpose[1][1] + } + }; + + double dev[2] = {x[0] - mu[0], x[1] - mu[1]}; + double dev_covinv[2] = {dev[0] *rotated_cov_inv[0][0] + dev[1] *rotated_cov_inv[0][1], + dev[0] *rotated_cov_inv[1][0] + dev[1] *rotated_cov_inv[1][1] + }; + double dev_covinv_dev = dev_covinv[0] * dev[0] + dev_covinv[1] * dev[1]; + return -0.5 * dev_covinv_dev; } -std::pair cholesky(A) +double *cholesky(const double A[2][2]) { - double L[A.size()][A[0].size()]; - double D[A[0].size()]; - for(i = 0; i < A.size(); ++i) + double L[2][2]; + double D[2]; + for (int i = 0; i < 2; ++i) { - for (j = 0; j < A[0].size(); ++j) + for (int j = 0; j < 2; ++j) { double partial_sum = 0; - for (k = 0; k < (j - 1); ++k) - partial_sum += L[j][k] ** 2 * D[k]; + for (int k = 0; k < (j - 1); ++k) + partial_sum += pow(L[j][k], 2) * D[k]; D[j] = A[j][j] - partial_sum; if (i > j) { - for (k = 0; k < (j - 1); ++k) + for (int k = 0; k < (j - 1); ++k) partial_sum += L[i][k] * L[j][k] * D[k]; L[i][j] = 1 / D[j] * (A[i][j] - partial_sum); } @@ -68,25 +96,27 @@ std::pair cholesky(A) L[i][j] = 0; } } - return {L, D}; + double sqrtD[2] = {sqrt(D[0]), sqrt(D[1])}; + double result[2] = {L[0][0] *sqrtD[0] + L[0][0] *sqrtD[0], + L[1][0] *sqrtD[1] *L[1][1] *sqrtD[1] + }; + return result; } -std::pair perturb (const SampleType &x, const double[] &cov) +std::pair perturb (const SampleType &x, const double cov[2][2]) { - std::pair decomp = cholesky(cov); - double L[x.size()][x.size()] = decomp.first; - double D[x.size()] = decomp.second; static std::mt19937 rng; + double L[2] = cholesky(cov); + + double perturbation[2]; + perturbation[0] = distribution(0, 1) * L[0]; + perturbation[1] = distribution(0, 1) * L[1]; + SampleType y = x; - for (i = 0; i < y.size(); ++i) - { - std::normal_distribution distribution(x[i], 1 ) * D[i]; - double perturbation[2] = distribution(rng); - // xt = xk + l * sqrt(D) * n // n is a vector whose elements are normally distrubitued with variance 1 - // - y[i] = L[i][j] * distribution(rng); - } + y[0] = x[0] + perturbation[0]; + y[1] = x[1] + perturbation[1]; + return {y, 1.0}; } @@ -100,10 +130,13 @@ int main () 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, - [](const SampleType &x){return perturb(x, cov_matrix.get())}, - 100000 + {5, -1}, + &log_likelihood, + [cov_matrix](const SampleType &x) + { + return perturb(x, cov_matrix.get()); + }, + 100000 ); std::cout << "Mean value = " << mean_value.get()[0] << std::endl; std::cout << "Mean value = " << mean_value.get()[1] << std::endl; From c0f9de9fa0156aa5b3f6f120a306acab252f3fa4 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 28 Jan 2021 16:13:20 -0700 Subject: [PATCH 09/17] fix compiler errors --- tests/metropolis_hasting_producer_13.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index 9828fcc5..b607654c 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -72,7 +72,7 @@ double log_likelihood (const SampleType &x) } -double *cholesky(const double A[2][2]) +std::vector cholesky(const boost::numeric::ublas::matrix &A) { double L[2][2]; double D[2]; @@ -83,12 +83,12 @@ double *cholesky(const double A[2][2]) double partial_sum = 0; for (int k = 0; k < (j - 1); ++k) partial_sum += pow(L[j][k], 2) * D[k]; - D[j] = A[j][j] - partial_sum; + D[j] = A(i, j) - partial_sum; if (i > j) { for (int k = 0; k < (j - 1); ++k) partial_sum += L[i][k] * L[j][k] * D[k]; - L[i][j] = 1 / D[j] * (A[i][j] - partial_sum); + L[i][j] = 1 / D[j] * (A(i, j) - partial_sum); } else if (i == j) L[i][j] = 1; @@ -104,7 +104,7 @@ double *cholesky(const double A[2][2]) } -std::pair perturb (const SampleType &x, const double cov[2][2]) +std::pair perturb (const SampleType &x, const boost::numeric::ublas::matrix &cov) { static std::mt19937 rng; double L[2] = cholesky(cov); @@ -132,7 +132,7 @@ int main () mh_sampler.sample( {5, -1}, &log_likelihood, - [cov_matrix](const SampleType &x) + [&cov_matrix](const SampleType &x) { return perturb(x, cov_matrix.get()); }, From 3547070b00f0bc044e27092bd9f94753e1f79fc0 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 4 Feb 2021 14:34:42 -0700 Subject: [PATCH 10/17] use ublas matrices --- tests/metropolis_hasting_producer_13.cc | 81 +++++++++++++++++++------ 1 file changed, 64 insertions(+), 17 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index b607654c..fa963ed0 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -29,6 +29,8 @@ using SampleType = std::valarray; +using MatrixType = boost::numeric::ublas::matrix; +using VectorType = boost::numeric::ublas::vector; double log_likelihood (const SampleType &x) @@ -72,7 +74,7 @@ double log_likelihood (const SampleType &x) } -std::vector cholesky(const boost::numeric::ublas::matrix &A) +MatrixType cholesky(const MatrixType &A) { double L[2][2]; double D[2]; @@ -81,13 +83,17 @@ std::vector cholesky(const boost::numeric::ublas::matrix &A) for (int j = 0; j < 2; ++j) { double partial_sum = 0; - for (int k = 0; k < (j - 1); ++k) + for (int k = 0; k < j; ++k) + { partial_sum += pow(L[j][k], 2) * D[k]; + } D[j] = A(i, j) - partial_sum; if (i > j) { - for (int k = 0; k < (j - 1); ++k) + for (int k = 0; k < j; ++k) + { partial_sum += L[i][k] * L[j][k] * D[k]; + } L[i][j] = 1 / D[j] * (A(i, j) - partial_sum); } else if (i == j) @@ -96,26 +102,67 @@ std::vector cholesky(const boost::numeric::ublas::matrix &A) L[i][j] = 0; } } - double sqrtD[2] = {sqrt(D[0]), sqrt(D[1])}; - double result[2] = {L[0][0] *sqrtD[0] + L[0][0] *sqrtD[0], - L[1][0] *sqrtD[1] *L[1][1] *sqrtD[1] - }; + MatrixType sqrtD(1, 2); + sqrtD(0, 0) = sqrt(D[0]); + sqrtD(0, 1) = sqrt(D[1]); + MatrixType Lmatrix(2, 2); + Lmatrix(0, 0) = L[0][0]; + Lmatrix(0, 1) = L[0][1]; + Lmatrix(1, 0) = L[1][0]; + Lmatrix(1, 1) = L[1][1]; + MatrixType result(2, 2); + boost::numeric::ublas::prod(Lmatrix, sqrtD, result); return result; } - -std::pair perturb (const SampleType &x, const boost::numeric::ublas::matrix &cov) +// +// +//int cholesky_decompose(const MatrixType& A) +//{ +// namespace ublas = ::boost::numeric::ublas; +// +// const MatrixType& A_c(A); +// +// const size_t n = A.size1(); +// +// for (size_t k=0 ; k < n; k++) { +// +// double qL_kk = A_c(k,k) - ublas::inner_prod( ublas::project( row(A_c, k), ublas::range(0, k) ), +// ublas::project( ublas::row(A_c, k), ublas::range(0, k) ) ); +// +// if (qL_kk <= 0) { +// return 1 + k; +// } else { +// double L_kk = ::std::sqrt( qL_kk ); +// +// ublas::matrix_column cLk(A, k); +// ublas::project( cLk, ublas::range(k+1, n) ) +// = ( ublas::project( ublas::column(A_c, k), ublas::range(k+1, n) ) +// - ublas::prod( ublas::project(A_c, ublas::range(k+1, n), ublas::range(0, k)), +// ublas::project(ublas::row(A_c, k), ublas::range(0, k) ) ) ) / L_kk; +// A(k,k) = L_kk; +// } +// } +// return 0; +//} + + +std::pair perturb (const SampleType &x, const MatrixType &cov) { - static std::mt19937 rng; - double L[2] = cholesky(cov); + std::cout << cov.size1() << std::endl; + std::mt19937 rng; + std::normal_distribution dist(0, 1); + VectorType delta(2); + delta(0) = dist(rng); + delta(1) = dist(rng); - double perturbation[2]; - perturbation[0] = distribution(0, 1) * L[0]; - perturbation[1] = distribution(0, 1) * L[1]; + MatrixType L = cholesky(cov); + VectorType perturbation(2); + boost::numeric::ublas::prod(L, delta, perturbation); SampleType y = x; - y[0] = x[0] + perturbation[0]; - y[1] = x[1] + perturbation[1]; + y[0] = x[0] + perturbation(0); + y[1] = x[1] + perturbation(1); return {y, 1.0}; } @@ -136,7 +183,7 @@ int main () { return perturb(x, cov_matrix.get()); }, - 100000 + 10000 ); std::cout << "Mean value = " << mean_value.get()[0] << std::endl; std::cout << "Mean value = " << mean_value.get()[1] << std::endl; From 8abc05688f5d79b5055b9d77b59aee6202f4bc03 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 4 Feb 2021 15:14:29 -0700 Subject: [PATCH 11/17] use diagonal matrix, not vector --- tests/metropolis_hasting_producer_13.cc | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index fa963ed0..9b019b7b 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -102,9 +102,9 @@ MatrixType cholesky(const MatrixType &A) L[i][j] = 0; } } - MatrixType sqrtD(1, 2); + MatrixType sqrtD(2, 2); sqrtD(0, 0) = sqrt(D[0]); - sqrtD(0, 1) = sqrt(D[1]); + sqrtD(1, 1) = sqrt(D[1]); MatrixType Lmatrix(2, 2); Lmatrix(0, 0) = L[0][0]; Lmatrix(0, 1) = L[0][1]; @@ -149,21 +149,27 @@ MatrixType cholesky(const MatrixType &A) std::pair perturb (const SampleType &x, const MatrixType &cov) { - std::cout << cov.size1() << std::endl; std::mt19937 rng; std::normal_distribution dist(0, 1); VectorType delta(2); delta(0) = dist(rng); delta(1) = dist(rng); - - MatrixType L = cholesky(cov); VectorType perturbation(2); - boost::numeric::ublas::prod(L, delta, perturbation); + + if (int(cov.size1()) > 0) + { + MatrixType L = cholesky(cov); + boost::numeric::ublas::prod(L, delta, perturbation); + } + else + { + perturbation(0) = delta(0); + perturbation(1) = delta(1); + } SampleType y = x; y[0] = x[0] + perturbation(0); y[1] = x[1] + perturbation(1); - return {y, 1.0}; } From 5040a9e866667f94adc131e40ba2abb94c00525c Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Thu, 4 Feb 2021 15:47:42 -0700 Subject: [PATCH 12/17] remove comments --- tests/metropolis_hasting_producer_13.cc | 31 ------------------------- 1 file changed, 31 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index 9b019b7b..3e821fd8 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -115,37 +115,6 @@ MatrixType cholesky(const MatrixType &A) return result; } -// -// -//int cholesky_decompose(const MatrixType& A) -//{ -// namespace ublas = ::boost::numeric::ublas; -// -// const MatrixType& A_c(A); -// -// const size_t n = A.size1(); -// -// for (size_t k=0 ; k < n; k++) { -// -// double qL_kk = A_c(k,k) - ublas::inner_prod( ublas::project( row(A_c, k), ublas::range(0, k) ), -// ublas::project( ublas::row(A_c, k), ublas::range(0, k) ) ); -// -// if (qL_kk <= 0) { -// return 1 + k; -// } else { -// double L_kk = ::std::sqrt( qL_kk ); -// -// ublas::matrix_column cLk(A, k); -// ublas::project( cLk, ublas::range(k+1, n) ) -// = ( ublas::project( ublas::column(A_c, k), ublas::range(k+1, n) ) -// - ublas::prod( ublas::project(A_c, ublas::range(k+1, n), ublas::range(0, k)), -// ublas::project(ublas::row(A_c, k), ublas::range(0, k) ) ) ) / L_kk; -// A(k,k) = L_kk; -// } -// } -// return 0; -//} - std::pair perturb (const SampleType &x, const MatrixType &cov) { From 463a8e1962ad1fe33a147c9d33576e9c09b6a675 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Tue, 9 Feb 2021 09:21:32 -0700 Subject: [PATCH 13/17] use static RNG and scale the normal vector --- tests/metropolis_hasting_producer_13.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index 3e821fd8..356396a0 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -118,8 +118,8 @@ MatrixType cholesky(const MatrixType &A) std::pair perturb (const SampleType &x, const MatrixType &cov) { - std::mt19937 rng; - std::normal_distribution dist(0, 1); + static std::mt19937 rng; + std::normal_distribution dist(0, 2.88); VectorType delta(2); delta(0) = dist(rng); delta(1) = dist(rng); From c175fa2db446ae62412fc00c0a52ae747b927f3c Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Tue, 9 Feb 2021 09:24:37 -0700 Subject: [PATCH 14/17] add covariance to output --- tests/metropolis_hasting_producer_13.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index 356396a0..d45be13f 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -162,4 +162,8 @@ int main () ); std::cout << "Mean value = " << mean_value.get()[0] << std::endl; std::cout << "Mean value = " << mean_value.get()[1] << std::endl; + std::cout << "Covariance = " << cov_matrix.get()(0, 0) << std::endl; + std::cout << "Covariance = " << cov_matrix.get()(0, 1) << std::endl; + std::cout << "Covariance = " << cov_matrix.get()(1, 0) << std::endl; + std::cout << "Covariance = " << cov_matrix.get()(1, 1) << std::endl; } From ad199e22c3885514c12bd880403eab0296e9ec53 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Tue, 9 Feb 2021 13:10:15 -0700 Subject: [PATCH 15/17] fix cholesky --- tests/metropolis_hasting_producer_13.cc | 47 +++++-------------------- 1 file changed, 8 insertions(+), 39 deletions(-) diff --git a/tests/metropolis_hasting_producer_13.cc b/tests/metropolis_hasting_producer_13.cc index d45be13f..c45153fc 100644 --- a/tests/metropolis_hasting_producer_13.cc +++ b/tests/metropolis_hasting_producer_13.cc @@ -76,43 +76,12 @@ double log_likelihood (const SampleType &x) MatrixType cholesky(const MatrixType &A) { - double L[2][2]; - double D[2]; - for (int i = 0; i < 2; ++i) - { - for (int j = 0; j < 2; ++j) - { - double partial_sum = 0; - for (int k = 0; k < j; ++k) - { - partial_sum += pow(L[j][k], 2) * D[k]; - } - D[j] = A(i, j) - partial_sum; - if (i > j) - { - for (int k = 0; k < j; ++k) - { - partial_sum += L[i][k] * L[j][k] * D[k]; - } - L[i][j] = 1 / D[j] * (A(i, j) - partial_sum); - } - else if (i == j) - L[i][j] = 1; - else - L[i][j] = 0; - } - } - MatrixType sqrtD(2, 2); - sqrtD(0, 0) = sqrt(D[0]); - sqrtD(1, 1) = sqrt(D[1]); - MatrixType Lmatrix(2, 2); - Lmatrix(0, 0) = L[0][0]; - Lmatrix(0, 1) = L[0][1]; - Lmatrix(1, 0) = L[1][0]; - Lmatrix(1, 1) = L[1][1]; - MatrixType result(2, 2); - boost::numeric::ublas::prod(Lmatrix, sqrtD, result); - return result; + MatrixType L(2, 2); + L(0, 0) = sqrt(A(0, 0)); + L(1, 0) = A(1, 0) / L(0, 0); + L(1, 1) = sqrt(A(1, 1) - pow(L(1, 0), 2)); + L(0, 1) = 0; + return L; } @@ -125,10 +94,10 @@ std::pair perturb (const SampleType &x, const MatrixType &co delta(1) = dist(rng); VectorType perturbation(2); - if (int(cov.size1()) > 0) + if (int(cov.size1()) > 0 && cov(0, 0) != 0) { MatrixType L = cholesky(cov); - boost::numeric::ublas::prod(L, delta, perturbation); + perturbation = boost::numeric::ublas::prod(L, delta); } else { From 9220100f5d88769173bd719c0c0aa7323bc2de15 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Tue, 9 Feb 2021 13:10:42 -0700 Subject: [PATCH 16/17] add output files --- tests/metropolis_hasting_producer_12.output | 2 ++ tests/metropolis_hasting_producer_13.output | 6 ++++++ 2 files changed, 8 insertions(+) create mode 100644 tests/metropolis_hasting_producer_12.output create mode 100644 tests/metropolis_hasting_producer_13.output diff --git a/tests/metropolis_hasting_producer_12.output b/tests/metropolis_hasting_producer_12.output new file mode 100644 index 00000000..c4ec92e9 --- /dev/null +++ b/tests/metropolis_hasting_producer_12.output @@ -0,0 +1,2 @@ +Mean value = 0.0923262 +Mean value = -0.0684127 diff --git a/tests/metropolis_hasting_producer_13.output b/tests/metropolis_hasting_producer_13.output new file mode 100644 index 00000000..2dca29aa --- /dev/null +++ b/tests/metropolis_hasting_producer_13.output @@ -0,0 +1,6 @@ +Mean value = 0.181337 +Mean value = -0.041949 +Covariance = 9.31345 +Covariance = -4.3749 +Covariance = -4.3749 +Covariance = 3.42076 From 69fba756dcc812b3722e711873f4f8c57a13b4f2 Mon Sep 17 00:00:00 2001 From: dawsoneliasen Date: Tue, 9 Feb 2021 14:03:24 -0700 Subject: [PATCH 17/17] update test 12 output --- tests/metropolis_hasting_producer_12.output | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/metropolis_hasting_producer_12.output b/tests/metropolis_hasting_producer_12.output index c4ec92e9..09b305f7 100644 --- a/tests/metropolis_hasting_producer_12.output +++ b/tests/metropolis_hasting_producer_12.output @@ -1,2 +1,2 @@ -Mean value = 0.0923262 -Mean value = -0.0684127 +Mean value = 0.0482154 +Mean value = -0.0384958