diff --git a/R/estimate.BSVARSIGN.R b/R/estimate.BSVARSIGN.R index ab19e9d..19405cd 100644 --- a/R/estimate.BSVARSIGN.R +++ b/R/estimate.BSVARSIGN.R @@ -96,16 +96,23 @@ estimate.BSVARSIGN = function(specification, S, thin = 1, show_progress = TRUE) # get the inputs to estimation # prior = specification$last_draw$prior$get_prior() - prior = specification$prior + prior = specification$prior$get_prior() starting_values = specification$starting_values$get_starting_values() identification = specification$identification$get_identification() max_tries = identification$max_tries data_matrices = specification$data_matrices$get_data_matrices() p = specification$p + + prior$B = t(prior$A) + prior$Ysoc = t(prior$Ysoc) + prior$Xsoc = t(prior$Xsoc) + prior$Ysur = t(prior$Ysur) + prior$Xsur = t(prior$Xsur) + Y = t(data_matrices$Y) + X = t(data_matrices$X) # estimation - qqq = .Call(`_bsvarSIGNs_bsvar_sign_cpp`, S, p, - data_matrices$Y, data_matrices$X, identification$VB, + qqq = .Call(`_bsvarSIGNs_bsvar_sign_cpp`, S, p, Y, X, identification$VB, identification$sign_irf, identification$sign_narrative, identification$sign_B, identification$zero_irf, prior, starting_values, show_progress, thin, max_tries) diff --git a/R/specify_bsvarSIGN.R b/R/specify_bsvarSIGN.R index 933f4b6..40ca716 100644 --- a/R/specify_bsvarSIGN.R +++ b/R/specify_bsvarSIGN.R @@ -139,23 +139,20 @@ specify_prior_bsvarSIGN = R6::R6Class( "PriorBSVARSIGN", public = list( + #' @field p a positive integer - the number of lags. + p = 1, + #' @field hyper a \code{(N+3)xS} matrix of hyper-parameters \eqn{\mu, \delta, \lambda, \psi}. hyper = matrix(), - #' @field B a \code{KxN} normal prior mean matrix for the autoregressive + #' @field A a \code{NxK} normal prior mean matrix for the autoregressive #' parameters. - B = matrix(), + A = matrix(), #' @field V a \code{KxK} matrix determining the normal prior column-specific #' covariance for the autoregressive parameters. V = matrix(), - #' @field Vp a \code{px1} vector of lag shrinkage level. - Vp = matrix(), - - #' @field Vd a \code{(d+1)x1} vector of (large) variances for constants. - Vd = matrix(), - #' @field S an \code{NxN} matrix determining the inverted-Wishart prior scale #' of error terms covariance matrix. S = matrix(), @@ -167,22 +164,22 @@ specify_prior_bsvarSIGN = R6::R6Class( #' @field data an \code{TxN} matrix of observations. data = matrix(), - #' @field Y an \code{TxN} matrix of dependent variables. + #' @field Y an \code{NxT} matrix of dependent variables. Y = matrix(), - #' @field X an \code{TxK} matrix of independent variables. + #' @field X an \code{KxT} matrix of independent variables. X = matrix(), #' @field Ysoc an \code{NxN} matrix with the sum-of-coefficients dummy observations. Ysoc = matrix(), - #' @field Xsoc an \code{NxK} matrix with the sum-of-coefficients dummy observations. + #' @field Xsoc an \code{KxN} matrix with the sum-of-coefficients dummy observations. Xsoc = matrix(), #' @field Ysur an \code{NxN} matrix with the single-unit-root dummy observations. Ysur = matrix(), - #' @field Xsur an \code{NxK} matrix with the single-unit-root dummy observations. + #' @field Xsur an \code{KxN} matrix with the single-unit-root dummy observations. Xsur = matrix(), #' @field mu.scale a positive scalar - the shape of the gamma prior for \eqn{\mu}. @@ -224,7 +221,7 @@ specify_prior_bsvarSIGN = R6::R6Class( #' prior = specify_prior_bsvarSIGN$new(oil, p = 1) #' prior$B # show autoregressive prior mean #' - initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, dim(data)[2])) { + initialize = function(data, p, exogenous = NULL, stationary = rep(FALSE, ncol(data))) { stopifnot("Argument p must be a positive integer number." = p > 0 & p %% 1 == 0) @@ -245,9 +242,7 @@ specify_prior_bsvarSIGN = R6::R6Class( B = matrix(0, K, N) B[1:N,] = diag(!stationary) - Vp = (1:p)^-2 - Vd = rep(100, 1 + d) - V = diag(c(kronecker(Vp, rep(1, N)), Vd)) + V = diag(c(kronecker((1:p)^-2, rep(1, N)), rep(100, 1 + d))) s2.ols = rep(NA, N) for (n in 1:N) { @@ -282,20 +277,18 @@ specify_prior_bsvarSIGN = R6::R6Class( # } # Xstar = cbind(Xstar, c(rep(0, N), 1), matrix(0, N + 1, d)) - self$Y = Y - self$X = X - self$Vp = Vp - self$Vd = Vd - + self$p = p self$hyper = hyper - self$B = B + self$A = t(B) self$V = V self$S = diag(N) self$nu = N + 2 - self$Ysoc = Ysoc - self$Xsoc = Xsoc - self$Ysur = Ysur - self$Xsur = Xsur + self$Y = t(Y) + self$X = t(X) + self$Ysoc = t(Ysoc) + self$Xsoc = t(Xsoc) + self$Ysur = t(Ysur) + self$Xsur = t(Xsur) self$mu.scale = scale self$mu.shape = shape self$delta.scale = scale @@ -316,11 +309,10 @@ specify_prior_bsvarSIGN = R6::R6Class( #' get_prior = function(){ list( + p = self$p, hyper = self$hyper, - B = self$B, + A = self$A, V = self$V, - Vp = self$Vp, - Vd = self$Vd, S = self$S, nu = self$nu, Ysoc = self$Ysoc, @@ -373,10 +365,16 @@ specify_prior_bsvarSIGN = R6::R6Class( init = narrow_hyper(model, hyper) prior = self$get_prior() + prior$B = t(prior$A) + prior$Ysoc = t(prior$Ysoc) + prior$Xsoc = t(prior$Xsoc) + prior$Ysur = t(prior$Ysur) + prior$Xsur = t(prior$Xsur) + result = stats::optim( init, \(x) -log_posterior_hyper(extend_hyper(hyper, model, matrix(x)), - model, self$Y, self$X, prior), + model, t(self$Y), t(self$X), prior), method = 'L-BFGS-B', lower = rep(0, length(init)), upper = init * 100, @@ -393,7 +391,8 @@ specify_prior_bsvarSIGN = R6::R6Class( variance = e$vectors %*% diag(as.vector(1 / abs(e$values))) %*% t(e$vectors) } - self$hyper = sample_hyper(S, burn_in, mode, model, self$Y, self$X, variance, prior) + self$hyper = sample_hyper(S, burn_in, mode, model, + t(self$Y), t(self$X), variance, prior) self$hyper = self$hyper[, -(1:burn_in)] } # END estimate_hyper @@ -716,15 +715,14 @@ specify_bsvarSIGN = R6::R6Class( B[lower.tri(B, diag = TRUE)] = TRUE self$data_matrices = bsvars::specify_data_matrices$new(data, p, exogenous) - self$data_matrices$Y = t(self$data_matrices$Y) - self$data_matrices$X = t(self$data_matrices$X) self$identification = specify_identification_bsvarSIGN$new(N, sign_irf, sign_narrative, sign_B, zero_irf, max_tries) - self$prior = specify_prior_bsvarSIGN$new(data, p, exogenous, stationary) + self$prior = specify_prior_bsvarSIGN$new(data, p, exogenous, + stationary) self$starting_values = bsvars::specify_starting_values_bsvar$new(N, self$p, d) }, # END initialize diff --git a/data/oil.rda b/data/oil.rda index ca5e5e8..cd1e14b 100644 Binary files a/data/oil.rda and b/data/oil.rda differ diff --git a/inst/include/bsvarSIGNs_RcppExports.h b/inst/include/bsvarSIGNs_RcppExports.h index 8580342..34ce1c7 100644 --- a/inst/include/bsvarSIGNs_RcppExports.h +++ b/inst/include/bsvarSIGNs_RcppExports.h @@ -88,7 +88,7 @@ namespace bsvarSIGNs { return Rcpp::as >(rcpp_result_gen); } - inline Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, const arma::mat& Y, const arma::mat& X, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const bool show_progress = true, const int thin = 100, const int& max_tries = 10000) { + inline Rcpp::List bsvar_sign_cpp(const int& S, const int& p, const arma::mat& Y, const arma::mat& X, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const bool show_progress = true, const int thin = 100, const int& max_tries = 10000) { typedef SEXP(*Ptr_bsvar_sign_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP); static Ptr_bsvar_sign_cpp p_bsvar_sign_cpp = NULL; if (p_bsvar_sign_cpp == NULL) { @@ -98,7 +98,7 @@ namespace bsvarSIGNs { RObject rcpp_result_gen; { RNGScope RCPP_rngScope_gen; - rcpp_result_gen = p_bsvar_sign_cpp(Shield(Rcpp::wrap(S)), Shield(Rcpp::wrap(lags)), Shield(Rcpp::wrap(Y)), Shield(Rcpp::wrap(X)), Shield(Rcpp::wrap(VB)), Shield(Rcpp::wrap(sign_irf)), Shield(Rcpp::wrap(sign_narrative)), Shield(Rcpp::wrap(sign_B)), Shield(Rcpp::wrap(Z)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(starting_values)), Shield(Rcpp::wrap(show_progress)), Shield(Rcpp::wrap(thin)), Shield(Rcpp::wrap(max_tries))); + rcpp_result_gen = p_bsvar_sign_cpp(Shield(Rcpp::wrap(S)), Shield(Rcpp::wrap(p)), Shield(Rcpp::wrap(Y)), Shield(Rcpp::wrap(X)), Shield(Rcpp::wrap(VB)), Shield(Rcpp::wrap(sign_irf)), Shield(Rcpp::wrap(sign_narrative)), Shield(Rcpp::wrap(sign_B)), Shield(Rcpp::wrap(Z)), Shield(Rcpp::wrap(prior)), Shield(Rcpp::wrap(starting_values)), Shield(Rcpp::wrap(show_progress)), Shield(Rcpp::wrap(thin)), Shield(Rcpp::wrap(max_tries))); } if (rcpp_result_gen.inherits("interrupted-error")) throw Rcpp::internal::InterruptedException(); diff --git a/inst/tinytest/test_specify.R b/inst/tinytest/test_specify.R index d2b740f..dfc7886 100644 --- a/inst/tinytest/test_specify.R +++ b/inst/tinytest/test_specify.R @@ -8,13 +8,13 @@ spec = specify_bsvarSIGN$new(oil) expect_identical(class(spec)[1], "BSVARSIGN") -expect_identical(dim(spec$data_matrices$Y)[1], - dim(spec$data_matrices$X)[1]) +expect_identical(dim(spec$data_matrices$Y)[2], + dim(spec$data_matrices$X)[2]) expect_identical(length(spec$data_matrices$get_data_matrices()), 2L) -expect_identical(dim(spec$prior$B)[1], dim(spec$data_matrices$X)[2]) +expect_identical(dim(spec$prior$A)[2], dim(spec$data_matrices$X)[1]) expect_identical(class(spec$prior$get_prior()), "list") @@ -51,13 +51,13 @@ spec = specify_bsvarSIGN$new(oil, expect_identical(class(spec)[1], "BSVARSIGN") -expect_identical(dim(spec$data_matrices$Y)[1], - dim(spec$data_matrices$X)[1]) +expect_identical(dim(spec$data_matrices$Y)[2], + dim(spec$data_matrices$X)[2]) expect_identical(length(spec$data_matrices$get_data_matrices()), 2L) -expect_identical(dim(spec$prior$B)[1], dim(spec$data_matrices$X)[2]) +expect_identical(dim(spec$prior$A)[2], dim(spec$data_matrices$X)[1]) expect_identical(class(spec$prior$get_prior()), "list") diff --git a/inst/varia/nicetry.R b/inst/varia/nicetry.R index 6b428e9..6344bae 100644 --- a/inst/varia/nicetry.R +++ b/inst/varia/nicetry.R @@ -4,7 +4,7 @@ name = sapply(1:7, \(i) data$ShortDescr[[i]][[1]]) data = data$y colnames(data) = name -data = optimism +# data = optimism spec = specify_bsvarSIGN$new(data, p = 4) @@ -12,9 +12,10 @@ start_time = Sys.time() spec$prior$estimate_hyper(S = 10000, burn_in = 5000, mu = TRUE, delta = TRUE, lambda = TRUE, psi = TRUE) -# post = estimate(spec, S = 1000) -# irf = compute_impulse_responses(post, horizon = 40) -# plot(irf, probability = 0.68) +post = estimate(spec, S = 100) +class(post) = "PosteriorBSVAR" +irf = compute_impulse_responses(post, horizon = 40) +plot(irf, probability = 0.68) end_time = Sys.time() end_time - start_time diff --git a/man/specify_prior_bsvarSIGN.Rd b/man/specify_prior_bsvarSIGN.Rd index 4d32ae4..f6ea882 100644 --- a/man/specify_prior_bsvarSIGN.Rd +++ b/man/specify_prior_bsvarSIGN.Rd @@ -45,18 +45,16 @@ prior$estimate_hyper(S = 5) \section{Public fields}{ \if{html}{\out{
}} \describe{ +\item{\code{p}}{a positive integer - the number of lags.} + \item{\code{hyper}}{a \code{(N+3)xS} matrix of hyper-parameters \eqn{\mu, \delta, \lambda, \psi}.} -\item{\code{B}}{a \code{KxN} normal prior mean matrix for the autoregressive +\item{\code{A}}{a \code{NxK} normal prior mean matrix for the autoregressive parameters.} \item{\code{V}}{a \code{KxK} matrix determining the normal prior column-specific covariance for the autoregressive parameters.} -\item{\code{Vp}}{a \code{px1} vector of lag shrinkage level.} - -\item{\code{Vd}}{a \code{(d+1)x1} vector of (large) variances for constants.} - \item{\code{S}}{an \code{NxN} matrix determining the inverted-Wishart prior scale of error terms covariance matrix.} @@ -65,17 +63,17 @@ inverted-Wishart prior for error terms covariance matrix.} \item{\code{data}}{an \code{TxN} matrix of observations.} -\item{\code{Y}}{an \code{TxN} matrix of dependent variables.} +\item{\code{Y}}{an \code{NxT} matrix of dependent variables.} -\item{\code{X}}{an \code{TxK} matrix of independent variables.} +\item{\code{X}}{an \code{KxT} matrix of independent variables.} \item{\code{Ysoc}}{an \code{NxN} matrix with the sum-of-coefficients dummy observations.} -\item{\code{Xsoc}}{an \code{NxK} matrix with the sum-of-coefficients dummy observations.} +\item{\code{Xsoc}}{an \code{KxN} matrix with the sum-of-coefficients dummy observations.} \item{\code{Ysur}}{an \code{NxN} matrix with the single-unit-root dummy observations.} -\item{\code{Xsur}}{an \code{NxK} matrix with the single-unit-root dummy observations.} +\item{\code{Xsur}}{an \code{KxN} matrix with the single-unit-root dummy observations.} \item{\code{mu.scale}}{a positive scalar - the shape of the gamma prior for \eqn{\mu}.} @@ -114,7 +112,7 @@ Create a new prior specification PriorBSVAR. data, p, exogenous = NULL, - stationary = rep(FALSE, dim(data)[2]) + stationary = rep(FALSE, ncol(data)) )}\if{html}{\out{
}} } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 85deade..7c44a68 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -127,12 +127,12 @@ RcppExport SEXP _bsvarSIGNs_bsvarSIGNs_ir(SEXP posterior_BSEXP, SEXP posterior_T return rcpp_result_gen; } // bsvar_sign_cpp -Rcpp::List bsvar_sign_cpp(const int& S, const int& lags, const arma::mat& Y, const arma::mat& X, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const bool show_progress, const int thin, const int& max_tries); -static SEXP _bsvarSIGNs_bsvar_sign_cpp_try(SEXP SSEXP, SEXP lagsSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) { +Rcpp::List bsvar_sign_cpp(const int& S, const int& p, const arma::mat& Y, const arma::mat& X, const arma::field& VB, const arma::cube& sign_irf, const arma::mat& sign_narrative, const arma::mat& sign_B, const arma::field& Z, const Rcpp::List& prior, const Rcpp::List& starting_values, const bool show_progress, const int thin, const int& max_tries); +static SEXP _bsvarSIGNs_bsvar_sign_cpp_try(SEXP SSEXP, SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::traits::input_parameter< const int& >::type S(SSEXP); - Rcpp::traits::input_parameter< const int& >::type lags(lagsSEXP); + Rcpp::traits::input_parameter< const int& >::type p(pSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); Rcpp::traits::input_parameter< const arma::field& >::type VB(VBSEXP); @@ -145,15 +145,15 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const bool >::type show_progress(show_progressSEXP); Rcpp::traits::input_parameter< const int >::type thin(thinSEXP); Rcpp::traits::input_parameter< const int& >::type max_tries(max_triesSEXP); - rcpp_result_gen = Rcpp::wrap(bsvar_sign_cpp(S, lags, Y, X, VB, sign_irf, sign_narrative, sign_B, Z, prior, starting_values, show_progress, thin, max_tries)); + rcpp_result_gen = Rcpp::wrap(bsvar_sign_cpp(S, p, Y, X, VB, sign_irf, sign_narrative, sign_B, Z, prior, starting_values, show_progress, thin, max_tries)); return rcpp_result_gen; END_RCPP_RETURN_ERROR } -RcppExport SEXP _bsvarSIGNs_bsvar_sign_cpp(SEXP SSEXP, SEXP lagsSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) { +RcppExport SEXP _bsvarSIGNs_bsvar_sign_cpp(SEXP SSEXP, SEXP pSEXP, SEXP YSEXP, SEXP XSEXP, SEXP VBSEXP, SEXP sign_irfSEXP, SEXP sign_narrativeSEXP, SEXP sign_BSEXP, SEXP ZSEXP, SEXP priorSEXP, SEXP starting_valuesSEXP, SEXP show_progressSEXP, SEXP thinSEXP, SEXP max_triesSEXP) { SEXP rcpp_result_gen; { Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = PROTECT(_bsvarSIGNs_bsvar_sign_cpp_try(SSEXP, lagsSEXP, YSEXP, XSEXP, VBSEXP, sign_irfSEXP, sign_narrativeSEXP, sign_BSEXP, ZSEXP, priorSEXP, starting_valuesSEXP, show_progressSEXP, thinSEXP, max_triesSEXP)); + rcpp_result_gen = PROTECT(_bsvarSIGNs_bsvar_sign_cpp_try(SSEXP, pSEXP, YSEXP, XSEXP, VBSEXP, sign_irfSEXP, sign_narrativeSEXP, sign_BSEXP, ZSEXP, priorSEXP, starting_valuesSEXP, show_progressSEXP, thinSEXP, max_triesSEXP)); } Rboolean rcpp_isInterrupt_gen = Rf_inherits(rcpp_result_gen, "interrupted-error"); if (rcpp_isInterrupt_gen) { diff --git a/src/bsvars_sign.cpp b/src/bsvars_sign.cpp index ead6630..e077562 100644 --- a/src/bsvars_sign.cpp +++ b/src/bsvars_sign.cpp @@ -17,7 +17,7 @@ using namespace arma; // [[Rcpp::export]] Rcpp::List bsvar_sign_cpp( const int& S, // number of draws from the posterior - const int& lags, // number of lags + const int& p, // number of lags const arma::mat& Y, // NxT dependent variables const arma::mat& X, // KxT dependent variables const arma::field& VB, // N-list @@ -56,7 +56,7 @@ Rcpp::List bsvar_sign_cpp( Rcout << " Press Esc to interrupt the computations" << endl; Rcout << "**************************************************|" << endl; } - Progress p(50, show_progress); + Progress bar(50, show_progress); const int T = Y.n_rows; const int N = Y.n_cols; @@ -79,6 +79,7 @@ Rcpp::List bsvar_sign_cpp( double w, mu, delta, lambda; vec hyper, psi; + vec prior_v = as(prior["V"]).diag(); mat B, Sigma, chol_Sigma, h_invp, Q, shocks; mat prior_V, prior_S, post_B, post_V, post_S; @@ -103,8 +104,14 @@ Rcpp::List bsvar_sign_cpp( psi = hyper.rows(3, N + 2); // update Minnesota prior - prior_V = diagmat(join_vert(lambda*lambda * kron(as(prior["Vp"]), 1 / psi), - as(prior["Vd"]))); + prior_v.rows(0, N * p - 1) = lambda * lambda * + prior_v.rows(0, N * p - 1) % + repmat(1 / psi, p, 1); + + prior_V = diagmat(lambda * lambda * prior_v % + join_vert(repmat(1 / psi, p, 1), + ones(K - N * p))); + prior_V = diagmat(prior_v); prior_S = diagmat(psi); // update dummy observation prior @@ -126,7 +133,7 @@ Rcpp::List bsvar_sign_cpp( B = rmatnorm_cpp(post_B, post_V, Sigma); h_invp = inv(trimatl(chol_Sigma)); // lower tri, h(Sigma) is upper tri - result = sample_Q(lags, Y, X, B, h_invp, chol_Sigma, prior, + result = sample_Q(p, Y, X, B, h_invp, chol_Sigma, prior, VB, sign_irf, sign_narrative, sign_B, Z, max_tries); Q = result(0); shocks = result(1); @@ -145,14 +152,14 @@ Rcpp::List bsvar_sign_cpp( if (s % 200 == 0) checkUserInterrupt(); // Increment progress bar - if (any(prog_rep_points == s)) p.increment(); + if (any(prog_rep_points == s)) bar.increment(); // if (omp_get_thread_num() == 0) { // // Check for user interrupts // if (s % 10 == 0) checkUserInterrupt(); // // // Increment progress bar - // if (any(prog_rep_points == s)) p.increment(); + // if (any(prog_rep_points == s)) bar.increment(); // } } // END s loop diff --git a/src/sample_hyper.cpp b/src/sample_hyper.cpp index 9e8ee58..70a0889 100644 --- a/src/sample_hyper.cpp +++ b/src/sample_hyper.cpp @@ -143,6 +143,7 @@ double log_ml_dummy( ) { int N = Y.n_cols; + int p = as(prior["p"]); double mu = hyper(0); double delta = hyper(1); double lambda = hyper(2); @@ -150,9 +151,9 @@ double log_ml_dummy( // update Minnesota prior mat prior_B = as(prior["B"]); - mat prior_V = diagmat(join_vert( - lambda*lambda * kron(as(prior["Vp"]), 1 / psi), - as(prior["Vd"]))); + vec v = as(prior["V"]).diag(); + v.rows(0, N*p - 1) = lambda * lambda * v.rows(0, N*p - 1) % repmat(1 / psi, p, 1); + mat prior_V = diagmat(v); mat prior_S = diagmat(psi); int prior_nu = as(prior["nu"]);