Skip to content

Commit

Permalink
corrected sampler for nu #6
Browse files Browse the repository at this point in the history
  • Loading branch information
donotdespair committed Jun 2, 2024
1 parent 44c5d30 commit cb5cbce
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 7 deletions.
2 changes: 1 addition & 1 deletion R/specify_bvarpanel.R
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ specify_starting_values_bvarPANEL = R6::R6Class(
self$A = matrix(stats::rnorm(K * N, sd = 0.001), K, N) + diag(K)[,1:N]
self$V = stats::rWishart(1, K + 1, diag(K))[,,1]
self$Sigma = stats::rWishart(1, N + 1, diag(N))[,,1]
self$nu = 30 #stats::rgamma(1, shape = 1, scale = 30)
self$nu = N + 1 + 0.1
self$m = stats::rnorm(1, sd = 0.001)
self$w = stats::rgamma(1, 1)
self$s = stats::rgamma(1, 1)
Expand Down
10 changes: 4 additions & 6 deletions src/sample_mniw.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,9 @@ double log_kernel_nu (

mat sum_aux_Sigma_c_inv(N, N);
for (int c = 0; c < C; c++) {
double ldSc = log_det_sympd(aux_Sigma_c_cpp.slice(c));
log_kernel_nu -= 0.5 * (aux_nu + N + K + 1) * ldSc;

sum_aux_Sigma_c_inv += aux_Sigma_c_inv.slice(c);
}
log_kernel_nu -= 0.5 * (aux_nu - N - 1) * trace(aux_Sigma * sum_aux_Sigma_c_inv);
Expand All @@ -144,11 +147,6 @@ double log_kernel_nu (
log_kernel_nu -= C * R::lgammafn(0.5 * (aux_nu + 1 - n));
} // EDN n loop

for (int c = 0; c < C; c++) {
double ldSc = log_det_sympd(aux_Sigma_c_cpp.slice(c));
log_kernel_nu -= 0.5 * (aux_nu + N + K + 1) * ldSc;
} // END c loop

return log_kernel_nu;
} // END log_kernel_nu

Expand Down Expand Up @@ -178,7 +176,7 @@ double sample_nu (
Cov_nu += R::psigamma( 0.5 * (aux_nu + 1 - n), 1);
} // END n loop
Cov_nu *= (C / 4);
Cov_nu -= (C * N * (N + 2)) * (2 * (aux_nu - N -1));
Cov_nu -= (C * N * aux_nu) * (2 * pow(aux_nu - N - 1, 2));
Cov_nu = sqrt(0.01 / Cov_nu);

// Metropolis-Hastings
Expand Down

0 comments on commit cb5cbce

Please sign in to comment.