Commit 291491d3 authored by chbauman's avatar chbauman
Browse files

uniform prior

parent f5960db7
......@@ -10,7 +10,7 @@
#include <omp.h>
#ifndef NDEBUG
#include "stefFenv.h"
//#include "stefFenv.h"
#endif
int main(){
......
......@@ -41,7 +41,7 @@ struct PotentialName<Potential::BUCKINGHAM> {
};
constexpr index_t POP_SIZE = 3000;
constexpr index_t POP_SIZE = 100000;
constexpr index_t N_DIM = (POTENTIAL == Potential::LENNARD_JONES? 2: 3) + 1;// + 1 for the unknown standard deviation
typedef Eigen::Matrix<numeric_t, N_DIM, 1> ThetaVector;
......@@ -61,3 +61,7 @@ constexpr numeric_t BETA = 0.2;
/// RNG
extern XoroshiroRandomNumberEngine GEN;
//#define UNIF_PRIOR
//#define USE_BRENTS_METHOD
constexpr numeric_t BOUND = 10e10;
......@@ -69,8 +69,18 @@ Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_exact_solution() {
constexpr numeric_t var_mean = (POTENTIAL == Potential::LENNARD_JONES? 4.823: 1.26075);
static std::normal_distribution<numeric_t> normal_prior(1.0, NORMAL_PRIOR_STDDEV); // 0.2 * 0.2
static std::gamma_distribution<numeric_t> var_dist(var_mean * var_mean, var_mean);
inline static PopMatrix prior_sample(const index_t n_samples){
static std::uniform_real_distribution<numeric_t> unif_dist(-1.,1.);
inline static PopMatrix prior_sample(const index_t n_samples){
PopMatrix prior(N_DIM, n_samples);
#ifdef UNIF_PRIOR
for(index_t j = 0; j < n_samples; ++j){
for(index_t i = 0; i < N_DIM - 1; ++i){
prior(i,j) = unif_dist(GEN) * BOUND;
}
prior(N_DIM - 1, j) = std::abs(unif_dist(GEN)) * BOUND;
}
#else
const auto exact_solution = get_exact_solution();
for(index_t i = 0; i < n_samples; ++i){
for (index_t j = 0; j < N_DIM-1; ++j) {
......@@ -78,14 +88,8 @@ inline static PopMatrix prior_sample(const index_t n_samples){
}
prior(N_DIM-1,i) = var_dist(GEN);
}
#endif
return prior;
/// For Buckingham:
for(index_t i = 0; i < N_DIM; ++i){
prior(0,i) = (0.2* normal_prior(GEN) + 1.0) * 1.45709e8;
prior(1,i) = (0.2* normal_prior(GEN) + 1.0) * 3.86461;
prior(2,i) = (0.2* normal_prior(GEN) + 1.0) * 468683;
prior(3,i) = var_dist(GEN);
}
}
/// Returns the value of the prior pdf at theta
......@@ -136,10 +140,19 @@ inline static PopMatrix mcmc(const ThetaVector & theta, const index_t rck, const
numeric_t particle_prop_likelihood;
particle_prop = sample_mvn(1, particle, sample_cov_mat);
particle_prop_likelihood = likelihood(particle_prop);
if(acc_dist(GEN) <= std::pow(particle_prop_likelihood / particle_likelihood, rho_curr) * prior(particle_prop) / prior(particle)){
#ifdef UNIF_PRIOR
if(particle_prop.cwiseAbs().maxCoeff() < BOUND &&
acc_dist(GEN) <= std::pow(particle_prop_likelihood / particle_likelihood, rho_curr)){
particle = particle_prop;
particle_likelihood = particle_prop_likelihood;
}
#else
if(acc_dist(GEN) <= std::pow(particle_prop_likelihood / particle_likelihood, rho_curr) *
prior(particle_prop) / prior(particle)){
particle = particle_prop;
particle_likelihood = particle_prop_likelihood;
}
#endif
new_thetas.col(i) = particle;
/*
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment