Skip to content

Commit

Permalink
correct sampler for nu closes #18
Browse files Browse the repository at this point in the history
  • Loading branch information
@donotdespair
donotdespair committed Jul 25, 2024
1 parent 2906ad3 commit 6afdbf8
Show file tree
Hide file tree
Showing 7 changed files with 137 additions and 80 deletions.
12 changes: 6 additions & 6 deletions R/RcppExports.R
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

bvarPANEL <- function(S, Y, X, prior, starting_values, thin, show_progress, rate_target_start_initial) {
.Call(`_bvarPANELs_bvarPANEL`, S, Y, X, prior, starting_values, thin, show_progress, rate_target_start_initial)
bvarPANEL <- function(S, Y, X, prior, starting_values, thin, show_progress, adptive_alpha_gamma) {
.Call(`_bvarPANELs_bvarPANEL`, S, Y, X, prior, starting_values, thin, show_progress, adptive_alpha_gamma)
}

Sigma2B_c <- function(posterior_Sigma_c, lower = TRUE) {
Expand Down Expand Up @@ -33,12 +33,12 @@ log_kernel_nu <- function(aux_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, p
.Call(`_bvarPANELs_log_kernel_nu`, aux_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior_lambda, C, N, K)
}

mcmc_accpetance_rate1 <- function(mcmc) {
.Call(`_bvarPANELs_mcmc_accpetance_rate1`, mcmc)
cov_nu <- function(aux_nu, C, N) {
.Call(`_bvarPANELs_cov_nu`, aux_nu, C, N)
}

sample_nu <- function(aux_nu, posterior_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior, iteration, scale, rate_target_start_initial) {
.Call(`_bvarPANELs_sample_nu`, aux_nu, posterior_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior, iteration, scale, rate_target_start_initial)
sample_nu <- function(aux_nu, adaptive_scale, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior, iteration, adptive_alpha_gamma) {
.Call(`_bvarPANELs_sample_nu`, aux_nu, adaptive_scale, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior, iteration, adptive_alpha_gamma)
}

sample_Sigma <- function(aux_Sigma_c_inv, aux_s, aux_nu, prior) {
Expand Down
2 changes: 1 addition & 1 deletion R/specify_bvarpanel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -441,7 +441,7 @@ specify_bvarPANEL = R6::R6Class(
self$data_matrices = specify_panel_data_matrices$new(data, self$p, exogenous)
self$prior = specify_prior_bvarPANEL$new(C, N, self$p, d)
self$starting_values = specify_starting_values_bvarPANEL$new(C, N, self$p, d)
self$adaptiveMH = c(0.6, 0.4, 10, 0.1)
self$adaptiveMH = c(0.44, 0.6)
}, # END initialize

#' @description
Expand Down
37 changes: 19 additions & 18 deletions src/RcppExports.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,8 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
#endif

// bvarPANEL
Rcpp::List bvarPANEL(const int& S, const Rcpp::List& Y, const Rcpp::List& X, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin, const bool show_progress, const arma::vec& rate_target_start_initial);
RcppExport SEXP _bvarPANELs_bvarPANEL(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP show_progressSEXP, SEXP rate_target_start_initialSEXP) {
Rcpp::List bvarPANEL(const int& S, const Rcpp::List& Y, const Rcpp::List& X, const Rcpp::List& prior, const Rcpp::List& starting_values, const int thin, const bool show_progress, const arma::vec& adptive_alpha_gamma);
RcppExport SEXP _bvarPANELs_bvarPANEL(SEXP SSEXP, SEXP YSEXP, SEXP XSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP thinSEXP, SEXP show_progressSEXP, SEXP adptive_alpha_gammaSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Expand All @@ -27,8 +27,8 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const Rcpp::List& >::type starting_values(starting_valuesSEXP);
Rcpp::traits::input_parameter< const int >::type thin(thinSEXP);
Rcpp::traits::input_parameter< const bool >::type show_progress(show_progressSEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type rate_target_start_initial(rate_target_start_initialSEXP);
rcpp_result_gen = Rcpp::wrap(bvarPANEL(S, Y, X, prior, starting_values, thin, show_progress, rate_target_start_initial));
Rcpp::traits::input_parameter< const arma::vec& >::type adptive_alpha_gamma(adptive_alpha_gammaSEXP);
rcpp_result_gen = Rcpp::wrap(bvarPANEL(S, Y, X, prior, starting_values, thin, show_progress, adptive_alpha_gamma));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -174,33 +174,34 @@ BEGIN_RCPP
return rcpp_result_gen;
END_RCPP
}
// mcmc_accpetance_rate1
double mcmc_accpetance_rate1(arma::vec& mcmc);
RcppExport SEXP _bvarPANELs_mcmc_accpetance_rate1(SEXP mcmcSEXP) {
// cov_nu
double cov_nu(const double& aux_nu, const int& C, const int& N);
RcppExport SEXP _bvarPANELs_cov_nu(SEXP aux_nuSEXP, SEXP CSEXP, SEXP NSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::vec& >::type mcmc(mcmcSEXP);
rcpp_result_gen = Rcpp::wrap(mcmc_accpetance_rate1(mcmc));
Rcpp::traits::input_parameter< const double& >::type aux_nu(aux_nuSEXP);
Rcpp::traits::input_parameter< const int& >::type C(CSEXP);
Rcpp::traits::input_parameter< const int& >::type N(NSEXP);
rcpp_result_gen = Rcpp::wrap(cov_nu(aux_nu, C, N));
return rcpp_result_gen;
END_RCPP
}
// sample_nu
double sample_nu(const double& aux_nu, const arma::vec& posterior_nu, const arma::cube& aux_Sigma_c_cpp, const arma::cube& aux_Sigma_c_inv, const arma::mat& aux_Sigma, const Rcpp::List& prior, const int& iteration, arma::vec& scale, const arma::vec& rate_target_start_initial);
RcppExport SEXP _bvarPANELs_sample_nu(SEXP aux_nuSEXP, SEXP posterior_nuSEXP, SEXP aux_Sigma_c_cppSEXP, SEXP aux_Sigma_c_invSEXP, SEXP aux_SigmaSEXP, SEXP priorSEXP, SEXP iterationSEXP, SEXP scaleSEXP, SEXP rate_target_start_initialSEXP) {
arma::vec sample_nu(double& aux_nu, double& adaptive_scale, const arma::cube& aux_Sigma_c_cpp, const arma::cube& aux_Sigma_c_inv, const arma::mat& aux_Sigma, const Rcpp::List& prior, const int& iteration, const arma::vec& adptive_alpha_gamma);
RcppExport SEXP _bvarPANELs_sample_nu(SEXP aux_nuSEXP, SEXP adaptive_scaleSEXP, SEXP aux_Sigma_c_cppSEXP, SEXP aux_Sigma_c_invSEXP, SEXP aux_SigmaSEXP, SEXP priorSEXP, SEXP iterationSEXP, SEXP adptive_alpha_gammaSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< const double& >::type aux_nu(aux_nuSEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type posterior_nu(posterior_nuSEXP);
Rcpp::traits::input_parameter< double& >::type aux_nu(aux_nuSEXP);
Rcpp::traits::input_parameter< double& >::type adaptive_scale(adaptive_scaleSEXP);
Rcpp::traits::input_parameter< const arma::cube& >::type aux_Sigma_c_cpp(aux_Sigma_c_cppSEXP);
Rcpp::traits::input_parameter< const arma::cube& >::type aux_Sigma_c_inv(aux_Sigma_c_invSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type aux_Sigma(aux_SigmaSEXP);
Rcpp::traits::input_parameter< const Rcpp::List& >::type prior(priorSEXP);
Rcpp::traits::input_parameter< const int& >::type iteration(iterationSEXP);
Rcpp::traits::input_parameter< arma::vec& >::type scale(scaleSEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type rate_target_start_initial(rate_target_start_initialSEXP);
rcpp_result_gen = Rcpp::wrap(sample_nu(aux_nu, posterior_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior, iteration, scale, rate_target_start_initial));
Rcpp::traits::input_parameter< const arma::vec& >::type adptive_alpha_gamma(adptive_alpha_gammaSEXP);
rcpp_result_gen = Rcpp::wrap(sample_nu(aux_nu, adaptive_scale, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior, iteration, adptive_alpha_gamma));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -277,8 +278,8 @@ static const R_CallMethodDef CallEntries[] = {
{"_bvarPANELs_sample_w", (DL_FUNC) &_bvarPANELs_sample_w, 2},
{"_bvarPANELs_sample_s", (DL_FUNC) &_bvarPANELs_sample_s, 5},
{"_bvarPANELs_log_kernel_nu", (DL_FUNC) &_bvarPANELs_log_kernel_nu, 8},
{"_bvarPANELs_mcmc_accpetance_rate1", (DL_FUNC) &_bvarPANELs_mcmc_accpetance_rate1, 1},
{"_bvarPANELs_sample_nu", (DL_FUNC) &_bvarPANELs_sample_nu, 9},
{"_bvarPANELs_cov_nu", (DL_FUNC) &_bvarPANELs_cov_nu, 3},
{"_bvarPANELs_sample_nu", (DL_FUNC) &_bvarPANELs_sample_nu, 8},
{"_bvarPANELs_sample_Sigma", (DL_FUNC) &_bvarPANELs_sample_Sigma, 4},
{"_bvarPANELs_sample_AV", (DL_FUNC) &_bvarPANELs_sample_AV, 6},
{"_bvarPANELs_sample_A_c_Sigma_c", (DL_FUNC) &_bvarPANELs_sample_A_c_Sigma_c, 6},
Expand Down
14 changes: 11 additions & 3 deletions src/bvarPANEL.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ Rcpp::List bvarPANEL(
const Rcpp::List& starting_values,
const int thin, // introduce thinning
const bool show_progress,
const arma::vec& rate_target_start_initial
const arma::vec& adptive_alpha_gamma // 2x1 vector with target acceptance rate and step size
) {

// Progress bar setup
Expand Down Expand Up @@ -82,6 +82,11 @@ Rcpp::List bvarPANEL(
vec scale(S);
int ss = 0;

// the initial value for the adaptive_scale is set to the negative inverse of
// Hessian for the posterior log_kenel for nu
double adaptive_scale = cov_nu(aux_nu, C, N);
vec aux_nu_tmp(2);

for (int s=0; s<S; s++) {
// Rcout << "Iteration: " << s << endl;

Expand All @@ -102,8 +107,11 @@ Rcpp::List bvarPANEL(

// sample aux_nu
// Rcout << " sample nu" << endl;
aux_nu = sample_nu( aux_nu, posterior_nu_S, aux_Sigma_c, aux_Sigma_c_inv, aux_Sigma, prior , s, scale, rate_target_start_initial);

// aux_nu = sample_nu( aux_nu, posterior_nu_S, aux_Sigma_c, aux_Sigma_c_inv, aux_Sigma, prior , s, scale, rate_target_start_initial);
aux_nu_tmp = sample_nu ( aux_nu, adaptive_scale, aux_Sigma_c, aux_Sigma_c_inv, aux_Sigma, prior, s, adptive_alpha_gamma );
aux_nu = aux_nu_tmp(0);
scale(s) = aux_nu_tmp(1);

// sample aux_Sigma
// Rcout << " sample Sigma" << endl;
aux_Sigma = sample_Sigma( aux_Sigma_c_inv, aux_s, aux_nu, prior );
Expand Down
2 changes: 1 addition & 1 deletion src/bvarPANEL.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ Rcpp::List bvarPANEL(
const Rcpp::List& starting_values,
const int thin, // introduce thinning
const bool show_progress,
const arma::vec& rate_target_start_initial
const arma::vec& adptive_alpha_gamma // 2x1 vector with target acceptance rate and step size
);


Expand Down
135 changes: 91 additions & 44 deletions src/sample_mniw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -153,38 +153,38 @@ double log_kernel_nu (




// [[Rcpp:interface(cpp)]]
// [[Rcpp::export]]
double mcmc_accpetance_rate1 (
arma::vec& mcmc
double cov_nu (
const double& aux_nu,
const int& C,
const int& N
) {
const int N = mcmc.n_elem;
double acc = 0;
for (int i=1; i<N; i++) {
if (mcmc(i) == mcmc(i-1)) {
acc++;
}
}
return 1 - (acc / (N - 1));
} // END mcmc_accpetance_rate1



double Cov_nu = 0;
for (int n = 1; n < N + 1; n++) {
Cov_nu += R::psigamma( 0.5 * (aux_nu + 1 - n), 1);
} // END n loop
Cov_nu *= (C / 4);
Cov_nu -= (C * N * aux_nu) * (2 * pow(aux_nu - N - 1, 2));
Cov_nu = sqrt(0.01 / Cov_nu);
return Cov_nu;
}



// [[Rcpp:interface(cpp)]]
// [[Rcpp::export]]
double sample_nu (
const double& aux_nu, // scalar
const arma::vec& posterior_nu, // sx1
arma::vec sample_nu (
double& aux_nu, // scalar
double& adaptive_scale,
const arma::cube& aux_Sigma_c_cpp, // NxNxC
const arma::cube& aux_Sigma_c_inv, // NxNxC
const arma::mat& aux_Sigma, // NxN
const Rcpp::List& prior,
const int& iteration, // MCMC iteration passed
arma::vec& scale, // (Sx1) adaptive scaling
const arma::vec& rate_target_start_initial
const arma::vec& adptive_alpha_gamma // 2x1 vector with target acceptance rate and step size
) {

double prior_lambda = as<double>(prior["lambda"]);
Expand All @@ -194,41 +194,88 @@ double sample_nu (
int N = aux_Sigma.n_rows;

// negative inverted Hessian of full conditional log-kernel
double Cov_nu = 0;
for (int n = 1; n < N + 1; n++) {
Cov_nu += R::psigamma( 0.5 * (aux_nu + 1 - n), 1);
} // END n loop
Cov_nu *= (C / 4);
Cov_nu -= (C * N * aux_nu) * (2 * pow(aux_nu - N - 1, 2));
Cov_nu = sqrt(0.01 / Cov_nu);

// Adaptive MH scaling
double scale_s = rate_target_start_initial(3);
if (iteration > rate_target_start_initial(2)) {
vec nu_to_s = posterior_nu.head(iteration - 1);
double alpha_s = mcmc_accpetance_rate1( nu_to_s );
scale_s = scale(iteration - 1) + pow(iteration, - rate_target_start_initial(0)) * (alpha_s - rate_target_start_initial(1));
}
scale(iteration) = scale_s;
double Cov_nu = cov_nu(aux_nu, C, N);

// Metropolis-Hastings
double aux_nu_star = RcppTN::rtn1( aux_nu, pow(scale_s, 2) * Cov_nu, N + 1, R_PosInf );
double aux_nu_star = RcppTN::rtn1( aux_nu, pow(adaptive_scale * Cov_nu, 0.5), N + 1, R_PosInf );
double lk_nu_star = log_kernel_nu ( aux_nu_star, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior_lambda, C, N, K );
double lk_nu_old = log_kernel_nu ( aux_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior_lambda, C, N, K );
double cgd_ratio = RcppTN::dtn1( aux_nu_star, aux_nu, pow(scale_s, 2) * Cov_nu, N + 1, R_PosInf ) /
RcppTN::dtn1( aux_nu, aux_nu_star, pow(scale_s, 2) * Cov_nu, N + 1, R_PosInf );
double cgd_ratio = RcppTN::dtn1( aux_nu_star, aux_nu, pow(adaptive_scale * Cov_nu, 0.5), N + 1, R_PosInf ) /
RcppTN::dtn1( aux_nu, aux_nu_star, pow(adaptive_scale * Cov_nu, 0.5), N + 1, R_PosInf );

double alpha = 1;
double kernel_ratio = exp(lk_nu_star - lk_nu_old) * cgd_ratio;
if ( kernel_ratio < 1 ) alpha = kernel_ratio;

double u = randu();
double out = 0;
if (u < exp(lk_nu_star - lk_nu_old) * cgd_ratio) {
out = aux_nu_star;
} else {
out = aux_nu;
} //
if ( u < alpha ) {
aux_nu = aux_nu_star;
}

if (iteration > 1) {
adaptive_scale = exp( log(adaptive_scale) + 0.5 * log( 1 + pow(iteration, - adptive_alpha_gamma(1)) * (alpha - adptive_alpha_gamma(0))) );
}

vec out = {aux_nu, adaptive_scale};
return out;
} // END sample_nu

// // [[Rcpp:interface(cpp)]]
// // [[Rcpp::export]]
// double sample_nu (
// const double& aux_nu, // scalar
// const arma::vec& posterior_nu, // sx1
// const arma::cube& aux_Sigma_c_cpp, // NxNxC
// const arma::cube& aux_Sigma_c_inv, // NxNxC
// const arma::mat& aux_Sigma, // NxN
// const Rcpp::List& prior,
// const int& iteration, // MCMC iteration passed
// arma::vec& scale, // (Sx1) adaptive scaling
// const arma::vec& rate_target_start_initial
// ) {
//
// double prior_lambda = as<double>(prior["lambda"]);
// mat prior_M = as<mat>(prior["M"]);
// int K = prior_M.n_rows;
// int C = aux_Sigma_c_cpp.n_slices;
// int N = aux_Sigma.n_rows;
//
// // negative inverted Hessian of full conditional log-kernel
// double Cov_nu = 0;
// for (int n = 1; n < N + 1; n++) {
// Cov_nu += R::psigamma( 0.5 * (aux_nu + 1 - n), 1);
// } // END n loop
// Cov_nu *= (C / 4);
// Cov_nu -= (C * N * aux_nu) * (2 * pow(aux_nu - N - 1, 2));
// Cov_nu = sqrt(0.01 / Cov_nu);
//
// // Adaptive MH scaling
// double scale_s = rate_target_start_initial(3);
// if (iteration > rate_target_start_initial(2)) {
// vec nu_to_s = posterior_nu.head(iteration - 1);
// double alpha_s = mcmc_accpetance_rate1( nu_to_s );
// scale_s = scale(iteration - 1) + pow(iteration, - rate_target_start_initial(0)) * (alpha_s - rate_target_start_initial(1));
// }
// scale(iteration) = scale_s;
//
// // Metropolis-Hastings
// double aux_nu_star = RcppTN::rtn1( aux_nu, pow(scale_s, 2) * Cov_nu, N + 1, R_PosInf );
// double lk_nu_star = log_kernel_nu ( aux_nu_star, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior_lambda, C, N, K );
// double lk_nu_old = log_kernel_nu ( aux_nu, aux_Sigma_c_cpp, aux_Sigma_c_inv, aux_Sigma, prior_lambda, C, N, K );
// double cgd_ratio = RcppTN::dtn1( aux_nu_star, aux_nu, pow(scale_s, 2) * Cov_nu, N + 1, R_PosInf ) /
// RcppTN::dtn1( aux_nu, aux_nu_star, pow(scale_s, 2) * Cov_nu, N + 1, R_PosInf );
//
// double u = randu();
// double out = 0;
// if (u < exp(lk_nu_star - lk_nu_old) * cgd_ratio) {
// out = aux_nu_star;
// } else {
// out = aux_nu;
// } //
//
// return out;
// } // END sample_nu



// [[Rcpp:interface(cpp)]]
Expand Down
15 changes: 8 additions & 7 deletions src/sample_mniw.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,21 +51,22 @@ double log_kernel_nu (
);


double mcmc_accpetance_rate1 (
arma::vec& mcmc
double cov_nu (
const double& aux_nu,
const int& C,
const int& N
);


double sample_nu (
const double& aux_nu, // scalar
const arma::vec& posterior_nu, // sx1
arma::vec sample_nu (
double& aux_nu, // scalar
double& adaptive_scale,
const arma::cube& aux_Sigma_c_cpp, // NxNxC
const arma::cube& aux_Sigma_c_inv, // NxNxC
const arma::mat& aux_Sigma, // NxN
const Rcpp::List& prior,
const int& iteration, // MCMC iteration passed
arma::vec& scale, // (Sx1) adaptive scaling
const arma::vec& rate_target_start_initial
const arma::vec& adptive_alpha_gamma // 2x1 vector with target acceptance rate and step size
);


Expand Down

0 comments on commit 6afdbf8

Please sign in to comment.