Commit 42a5e6fb authored by stefanow's avatar stefanow
Browse files

Faster sampling

parent a4ad14f4
......@@ -2,22 +2,22 @@ static std::normal_distribution<> dist_mvn;
struct normal_random_variable
{
normal_random_variable(Eigen::MatrixXd const& covar)
: normal_random_variable(Eigen::VectorXd::Zero(covar.rows()), covar)
{}
normal_random_variable(Eigen::VectorXd const& mean, Eigen::MatrixXd const& covar)
: mean(mean)
normal_random_variable(Eigen::MatrixXd const& covar)
{
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(covar);
transform = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
}
Eigen::VectorXd mean;
Eigen::MatrixXd transform;
Eigen::VectorXd operator()() const
{
return mean + transform * Eigen::VectorXd{ mean.size() }.unaryExpr([&](auto x) { return dist_mvn(GEN); });
const int numberOfEntries = transform.cols();
Eigen::VectorXd randomStandardNormal(numberOfEntries);
for (int i = 0; i < numberOfEntries; ++i) {
randomStandardNormal(i) = dist_mvn(GEN);
}
return transform * randomStandardNormal;
}
};
......@@ -63,7 +63,7 @@ Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_exact_solution() {
return exact_solution;
}
/// Returns sample from prior f(\Theta)
constexpr numeric_t var_mean = (POTENTIAL == Potential::LENNARD_JONES? 4.823: 1.26075);
......@@ -106,11 +106,6 @@ inline static numeric_t prior(const ThetaVector & theta){
return prior_pdf;
}
/// Returns 'n_samples' samples from a MVN distribution
inline static PopMatrix sample_mvn(const index_t n_samples, const ThetaVector &meanVec, const CovMat &covMat){
normal_random_variable sample { meanVec, covMat };
return sample();
};
/// Does the resampling
inline static void resample(const WeightVector & weights, CountVector & resample_count){
......@@ -168,7 +163,7 @@ inline static PopMatrix mcmc(const normal_random_variable & sample, const ThetaV
new_likelihoods(ct + i) = particle_likelihood;
//std::cout << particle_likelihood << "<-part. lh\n";
#endif
/*
/*
std::cout << "Cov mat " << sample_cov_mat << "\n"
<< "cur part " << particle << "\n"
<< "prop part " << particle_prop << "\n"
......
Markdown is supported
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