# [Rcpp-devel] rwishart Rcpp implementation

Yue Li gorillayue at gmail.com
Thu Jun 11 21:59:28 CEST 2015

```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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150611/8ce2a5a6/attachment.html>
```