[Rcpp-devel] rwishart Rcpp implementation

[Neal Fultz]

I had ported the bayesm version a while a back, here it is:

List rwishart(int nu,  NumericMatrix const& V){// function to draw
from Wishart (nu,V) and IW// // W ~ W(nu,V) E[W]=nuV// // WI=W^-1
E[WI]=V^-1/(nu-m-1)
  RNGScope rngscope;  int m = V.nrow();    mat Vm(V.begin(), m, m,
false);    // Can't vectorise because Rcpp-sugar rchisq doesnt
vectorise df argument  // T has sqrt chisqs on diagonal and normals
below diagonal and garbage above diagonal  mat T(m,m);    for(int i =
0; i < m; i++) {    T(i,i) = sqrt(R::rchisq(nu-i));  }  for(int j = 0;
j < m; j++) {    for(int i = j+1; i < m; i++) {          T(i,j) =
R::norm_rand();  }}    // explicitly declare T as triangular  // top
triangular is NaN ###  mat C = trans(trimatl(T)) * chol(Vm);  mat CI =
C.i();// C is the upper triangular root of Wishart // therefore, W=C'C
 this is the LU decomposition // Inv(W) = CICI'    this is the UL
decomp *not LU*!// W is Wishart draw, IW is W^-1  return List::create(
   _["W"] = trans(C) * C,    _["IW"] = CI * trans(CI),    _["C"] = C,
  _["CI"] = CI    );}


On Thu, Jun 11, 2015 at 12:59 PM, Yue Li <gorillayue at gmail.com> wrote:

> 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
>
> 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 at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150611/ce0f7365/attachment.html>