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 <gorilla...@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@lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
_______________________________________________
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
