[Rcpp-devel] rwishart Rcpp implementation

[Yue Li]

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>