Skip to content

Commit 0a98409

Browse files
zhizuiozhizuio
committed
disable births/deaths for RW-MH
1 parent 7f5144c commit 0a98409

File tree

2 files changed

+34
-19
lines changed

2 files changed

+34
-19
lines changed

‎src/BVS.cpp‎

Lines changed: 32 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -396,14 +396,15 @@ void BVS_Sampler::sampleGamma(
396396
// iterations for the leapfrog steps
397397
unsigned int I = 3;
398398
// add ΔH contributing to the Hastings correction
399-
logProposalRatio += HamiltonianBetas(I, proposedBeta, updateIdx0, componentUpdateIdx,
399+
/*logProposalRatio +=*/ (void)HamiltonianBetas(I, proposedBeta, betas_, updateIdx0, componentUpdateIdx,
400400
datTheta, datProportion, weibullS, weibullLambda, kappa_, tauSq_[componentUpdateIdx], sigmaMH_beta, dataclass);
401+
//volume-preserving and reversible HMC needs no proposal term when you do its own accept/reject at the kernel level.
401402
}
402403
}
403404

404405
} else {//if( arma::any(gammas_(updateIdx,singleIdx_k)) ) {
405406
// (symmetric) random-walk Metropolis with optimal standard deviation O(d^{-1/2}, theoretically 2.38*d^{-1/2})
406-
407+
/*
407408
// Symmetric RW-MH with birth handling
408409
arma::uvec birth = setdiff_preserve_order(updateIdx0, updateIdx0_rev); // newly active
409410
arma::uvec death = setdiff_preserve_order(updateIdx0_rev, updateIdx0); // becoming inactive
@@ -423,10 +424,13 @@ void BVS_Sampler::sampleGamma(
423424
arma::vec bcurr = betas_(1 + death, singleIdx_k); // current coefficients
424425
logProposalRatio += logPDFNormal(bcurr, s_birth_rev * s_birth_rev); // +log q(current|proposed)
425426
}
427+
unsigned int J_stay = stay.n_elem;
428+
*/
426429

427430
// Symmetric RW on stay only
428-
unsigned int J_stay = stay.n_elem;
431+
unsigned int J_stay = updateIdx0.n_elem;
429432
if (J_stay > 0) {
433+
arma::uvec stay = updateIdx0;
430434
double s_rw = sigmaMH_beta * 2.38 / std::sqrt(J_stay);
431435
arma::vec u = Rcpp::rnorm(J_stay, 0.0, s_rw);
432436
proposedBeta(1 + stay, singleIdx_k) += u;
@@ -440,8 +444,8 @@ void BVS_Sampler::sampleGamma(
440444
double logPriorBetaRatio = 0.;
441445
double logLikelihoodRatio = 0.;
442446
arma::vec proposedLikelihood = log_likelihood_;
443-
if(!used_hmc) // HMC already accounted in ΔH
444-
{
447+
// if(!used_hmc) // HMC already accounted in ΔH
448+
// {
445449
// prior ratio of beta
446450
logPriorBetaRatio += logPDFNormal(proposedBeta(1+updateIdx0,singleIdx_k), tauSq_[componentUpdateIdx]);
447451
logPriorBetaRatio -= logPDFNormal(betas_(1+updateIdx0_rev,singleIdx_k), tauSq_[componentUpdateIdx]);
@@ -450,7 +454,7 @@ void BVS_Sampler::sampleGamma(
450454
// loglikelihood( xi_, zetas_, betas_, kappa_, proportion_model, dataclass, log_likelihood_ ); // TODO: this can be removed
451455
loglikelihood( xi_, zetas_, proposedBeta, kappa_, proportion_model, dataclass, proposedLikelihood );
452456
logLikelihoodRatio = arma::sum(proposedLikelihood - log_likelihood_);
453-
}
457+
// }
454458

455459
// Here we need always compute the proposal and original ratios, in particular the likelihood, since betas are updated
456460
//logProposalGammaRatio = arma::accu(proposedGammaPrior - logP_gamma);
@@ -485,11 +489,11 @@ void BVS_Sampler::sampleGamma(
485489
}
486490
betas_ = proposedBeta;
487491
// Update log_likelihood_ at the accepted state
488-
if (used_hmc) {
489-
loglikelihood(xi_, zetas_, betas_, kappa_, proportion_model, dataclass, log_likelihood_);
490-
} else {
492+
// if (used_hmc) {
493+
// loglikelihood(xi_, zetas_, betas_, kappa_, proportion_model, dataclass, log_likelihood_);
494+
// } else {
491495
log_likelihood_ = proposedLikelihood;
492-
}
496+
// }
493497
// logPosteriorBeta = logPosteriorBeta_proposal;
494498

495499
++gamma_acc_count_;
@@ -750,13 +754,15 @@ void BVS_Sampler::sampleEta(
750754
// iterations for the leapfrog steps
751755
unsigned int I = 3;
752756
// add ΔH contributing to the Hastings correction
753-
logProposalRatio += HamiltonianZetas(I, proposedZeta, updateIdx0, componentUpdateIdx,
757+
/*logProposalRatio +=*/ (void)HamiltonianZetas(I, proposedZeta, zetas_, updateIdx0, componentUpdateIdx,
754758
datTheta, weibullS, weibullLambda, kappa_, wSq_[componentUpdateIdx], sigmaMH_zeta, dataclass);
759+
// volume-preserving and reversible HMC needs no proposal term when you do its own accept/reject at the kernel level.
755760
}
756761
}
757762

758763
} else {//if( arma::any(etas_(updateIdx,singleIdx_k)) ) {
759764
// (symmetric) random-walk Metropolis with optimal standard deviation O(d^{-1/2}, theoretically 2.38*d^{-1/2})
765+
/*
760766
// Symmetric RW-MH with birth handling
761767
arma::uvec birth = setdiff_preserve_order(updateIdx0, updateIdx0_rev);
762768
arma::uvec death = setdiff_preserve_order(updateIdx0_rev, updateIdx0);
@@ -776,10 +782,13 @@ void BVS_Sampler::sampleEta(
776782
arma::vec zcurr = zetas_(1 + death, singleIdx_k);
777783
logProposalRatio += logPDFNormal(zcurr, s_birth_rev * s_birth_rev);
778784
}
785+
unsigned int J_stay = stay.n_elem;
786+
*/
779787

780788
// Symmetric RW on stay
781-
unsigned int J_stay = stay.n_elem;
789+
unsigned int J_stay = updateIdx0.n_elem;
782790
if (J_stay > 0) {
791+
arma::uvec stay = updateIdx0;
783792
double s_rw = sigmaMH_zeta * 2.38 / std::sqrt(J_stay);
784793
arma::vec u = Rcpp::rnorm(J_stay, 0.0, s_rw);
785794
proposedZeta(1 + stay, singleIdx_k) += u;
@@ -792,8 +801,8 @@ void BVS_Sampler::sampleEta(
792801
double logPriorZetaRatio = 0.;
793802
double logLikelihoodRatio = 0.;
794803
arma::vec proposedLikelihood = log_likelihood_;
795-
if(!used_hmc) // HMC already accounted in ΔH
796-
{
804+
// if(!used_hmc) // HMC already accounted in ΔH
805+
// {
797806
logPriorZetaRatio += logPDFNormal(proposedZeta(1+updateIdx0,singleIdx_k), wSq_[componentUpdateIdx]);
798807
logPriorZetaRatio -= logPDFNormal(zetas_(1+updateIdx0_rev,singleIdx_k), wSq_[componentUpdateIdx]);
799808

@@ -802,7 +811,7 @@ void BVS_Sampler::sampleEta(
802811
// loglikelihood( xi_, zetas_, betas_, kappa_, true, dataclass, log_likelihood_ ); // TODO: this can be removed
803812
loglikelihood( xi_, proposedZeta, betas_, kappa_, true, dataclass, proposedLikelihood );
804813
logLikelihoodRatio = arma::sum(proposedLikelihood - log_likelihood_);
805-
}
814+
// }
806815

807816

808817
// Here we need always compute the proposal and original ratios, in particular the likelihood, since betas are updated
@@ -821,11 +830,11 @@ void BVS_Sampler::sampleEta(
821830

822831
}
823832
zetas_ = proposedZeta;
824-
if (used_hmc) {
825-
loglikelihood(xi_, zetas_, betas_, kappa_, true, dataclass, log_likelihood_);
826-
} else {
833+
// if (used_hmc) {
834+
// loglikelihood(xi_, zetas_, betas_, kappa_, true, dataclass, log_likelihood_);
835+
// } else {
827836
log_likelihood_ = proposedLikelihood;
828-
}
837+
// }
829838
// logPosteriorZeta = logPosteriorZeta_proposal;
830839

831840
++eta_acc_count_;
@@ -1229,6 +1238,7 @@ double BVS_Sampler::logPzetaK(
12291238
double BVS_Sampler::HamiltonianBetas(
12301239
unsigned int I,
12311240
arma::mat& proposedBeta,
1241+
const arma::mat& betas_,
12321242
const arma::uvec& updateIdx0,
12331243
unsigned int componentUpdateIdx,
12341244

@@ -1270,6 +1280,7 @@ double BVS_Sampler::HamiltonianBetas(
12701280
U = - logPbetaK(componentUpdateIdx, proposedBeta, tauSqK, kappa_, datTheta, datProportion, dataclass);
12711281
H -= U + 0.5 * arma::as_scalar(r.t() * r); // HMC target ratio (difference between intial and final energy)
12721282

1283+
if( std::log(R::runif(0,1)) >= H ) proposedBeta = betas_;
12731284
/*
12741285
if( std::log(R::runif(0,1)) < H )
12751286
{
@@ -1371,6 +1382,7 @@ arma::vec BVS_Sampler::gradientLogPbetas(
13711382
double BVS_Sampler::HamiltonianZetas(
13721383
unsigned int I,
13731384
arma::mat& proposedZeta,
1385+
const arma::mat& zetas_,
13741386
const arma::uvec& updateIdx0,
13751387
unsigned int componentUpdateIdx,
13761388

@@ -1411,6 +1423,7 @@ double BVS_Sampler::HamiltonianZetas(
14111423
U = -logPzetaK(componentUpdateIdx, proposedZeta, wSqK, kappa_, datTheta, weibullS, weibullLambda, dataclass);
14121424
H -= U + 0.5 * arma::as_scalar(r.t() * r); // HMC target ratio (difference between intial and final energy)
14131425

1426+
if( std::log(R::runif(0,1)) >= H ) proposedZeta = zetas_;
14141427
/*
14151428
if( std::log(R::runif(0,1)) < H )
14161429
{

‎src/BVS.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -215,6 +215,7 @@ class BVS_Sampler
215215
static double HamiltonianBetas(
216216
unsigned int I,
217217
arma::mat& proposedBeta,
218+
const arma::mat& betas_,
218219
const arma::uvec& updateIdx0,
219220
unsigned int componentUpdateIdx,
220221

@@ -243,6 +244,7 @@ class BVS_Sampler
243244
static double HamiltonianZetas(
244245
unsigned int I,
245246
arma::mat& proposedZeta,
247+
const arma::mat& zetas_,
246248
const arma::uvec& updateIdx0,
247249
unsigned int componentUpdateIdx,
248250

0 commit comments

Comments
 (0)