<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Dear List,<div class=""><br class=""></div><div class="">I wonder if anyone could share their Rcpp code for an equivalence of rWishart. I’ve come across a similar question in the forum:</div><div class=""><br class=""></div><div class=""><a href="http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp" class="">http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp</a></div><div class=""><br class=""></div><div class="">Based on Dirk’s suggestion, I have my following attempt:</div><div class=""><br class=""></div><div class=""><div class="">// [[Rcpp::export]]</div><div class="">mat mvrnormArma(int n, vec mu, mat sigma) {</div><div class=""> </div><div class=""> int ncols = sigma.n_cols;</div><div class=""> </div><div class=""> mat Y = randn(n, ncols);</div><div class=""> </div><div class=""> return repmat(mu, 1, n).t() + Y * chol(sigma);</div><div class="">}</div><div class=""><br class=""></div><div class=""><br class=""></div><div class="">// [[Rcpp::export]]</div><div class="">mat rWishartCpp(mat sigma, int n) {</div><div class=""> </div><div class=""> mat X = mvrnormArma(1, zeros<vec>(sigma.n_rows), sigma);</div><div class=""> </div><div class=""> return X.t() * X;</div><div class="">}</div></div><div class=""><br class=""></div><div class="">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:</div><div class=""><br class=""></div><div class="">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 …</div><div class=""><br class=""></div><div class=""><div class="">function baysm_rwishart(nu, V)</div><div class="">{</div><div class=""> m = nrow(V)</div><div class=""> </div><div class=""> df = (nu + nu - m + 1) - (nu - m + 1):nu</div><div class=""> </div><div class=""> if (m > 1) {</div><div class=""> T = diag(sqrt(rchisq(c(rep(1, m)), df)))</div><div class=""> T[lower.tri(T)] = rnorm((m * (m + 1)/2 - m))</div><div class=""> } </div><div class=""> else {</div><div class=""> T = sqrt(rchisq(1, df))</div><div class=""> }</div><div class=""> </div><div class=""> U = chol(V)</div><div class=""> </div><div class=""> C = t(T) %*% U</div><div class=""> </div><div class=""> CI = backsolve(C, diag(m))</div><div class=""> </div><div class=""> return(list(W = crossprod(C), </div><div class=""> <span class="Apple-tab-span" style="white-space:pre"> </span>IW = crossprod(t(CI)), C = C, </div><div class=""> CI = CI))</div><div class="">}</div></div><div class=""><br class=""></div><div class="">Yue</div><div class=""><br class=""></div><div class=""><br class=""></div></body></html>