[Rcpp-devel] rwishart Rcpp implementation

[Yue Li]

Many thanks Neal!

Yue

> On Jun 11, 2015, at 5:07 PM, Neal Fultz <nfultz at gmail.com> wrote:
> 
> 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 <mailto: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 <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 <mailto:Rcpp-devel at lists.r-forge.r-project.org>
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel <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/2ef536da/attachment-0001.html>