Commit 0d32d8ec authored by chbauman's avatar chbauman
Browse files

new sampling

parent f2e650e5
......@@ -99,7 +99,7 @@ int main(){
#endif
//std::cout << resample_count.transpose() << "<- res count\n";
//std::cout << indices_for_mcmc.transpose() << "<- indices\n";
/// Compute sample mean and covariance
sample_mean = thetas * weights;
......@@ -122,7 +122,7 @@ int main(){
const index_t rck = indices_for_mcmc(i, 0);
thetas_new.middleCols(ct, rck) = mcmc(sample, thetas.col(k), rck,
sample_mean, sample_cov_mat, rho_curr, likelihoods, likelihoods_old, likelihood, k, ct);
}
}
#else
index_t ct = 0;
#pragma omp parallel shared(sample_mean, sample_cov_mat, rho_curr, thetas)
......
......@@ -118,12 +118,14 @@ inline static PartPopMatrix mcmc(const normal_random_variable & sample, const Th
particle = theta;
numeric_t particle_likelihood = likelihoods(k);
#ifdef BURNIN
PartPopMatrix all_samples = sample(rck + BURNIN);
for(index_t i = 0; i < rck + BURNIN; ++i){
#else
PartPopMatrix all_samples = sample(rck);
for(index_t i = 0; i < rck; ++i){
#endif
numeric_t particle_prop_likelihood;
particle_prop = particle + sample();
particle_prop = particle + all_samples.col(i);
#ifdef UNIF_PRIOR
if(particle_prop.cwiseAbs().maxCoeff() <= BOUND && particle_prop.minCoeff() >= 0){
particle_prop_likelihood = likelihood(particle_prop);
......
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