Many thanks Neal!
Yue
> On Jun 11, 2015, at 5:07 PM, Neal Fultz <nfu...@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 <gorilla...@gmail.com
> <mailto: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
>
> <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
> <mailto:Rcpp-devel@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>
>
_______________________________________________
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
