Dear List, I wonder if anyone could share their Rcpp code for an equivalence of rWishart. I’ve come across a similar question in the forum:
http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp <http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp> Based on Dirk’s suggestion, I have my following attempt: // [[Rcpp::export]] mat mvrnormArma(int n, vec mu, mat sigma) { int ncols = sigma.n_cols; mat Y = randn(n, ncols); return repmat(mu, 1, n).t() + Y * chol(sigma); } // [[Rcpp::export]] mat rWishartCpp(mat sigma, int n) { mat X = mvrnormArma(1, zeros<vec>(sigma.n_rows), sigma); return X.t() * X; } But this implementation does not take into account the degree of freedom. In fact, I’m very fuzzy on how to incorporate df into the sampling process: I also checked the R code from baysm, namely rwishart, which samples from rchisq instead (code pasted below). But I’m not sure how to relate the simple naive code I have above with the baysm code below … function baysm_rwishart(nu, V) { m = nrow(V) df = (nu + nu - m + 1) - (nu - m + 1):nu if (m > 1) { T = diag(sqrt(rchisq(c(rep(1, m)), df))) T[lower.tri(T)] = rnorm((m * (m + 1)/2 - m)) } else { T = sqrt(rchisq(1, df)) } U = chol(V) C = t(T) %*% U CI = backsolve(C, diag(m)) return(list(W = crossprod(C), IW = crossprod(t(CI)), C = C, CI = CI)) } Yue
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
