<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Many thanks Neal!<div class=""><br class=""></div><div class="">Yue</div><div class=""><br class=""><div><blockquote type="cite" class=""><div class="">On Jun 11, 2015, at 5:07 PM, Neal Fultz <<a href="mailto:nfultz@gmail.com" class="">nfultz@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class="">I had ported the bayesm version a while a back, here it is:<div class=""><br class=""></div><div class=""><pre style="margin-top:0px;margin-bottom:0px;padding:0px;font-family:Consolas,Menlo,'Liberation Mono',Courier,monospace;font-size:12px;line-height:1.4;color:rgb(51,51,51)" class=""><span class="">List</span> <span class="" style="color:rgb(153,0,0);font-weight:bold">rwishart</span><span class="">(</span><span class="" style="color:rgb(68,85,136);font-weight:bold">int</span> <span class="">nu</span><span class="">,</span> <span class="">NumericMatrix</span> <span class="" style="font-weight:bold">const</span><span class="" style="font-weight:bold">&</span> <span class="">V</span><span class="">){</span>
<a name="cl-38" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// function to draw from Wishart (nu,V) and IW</span>
<a name="cl-39" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// </span>
<a name="cl-40" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// W ~ W(nu,V) E[W]=nuV</span>
<a name="cl-41" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// </span>
<a name="cl-42" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// WI=W^-1 E[WI]=V^-1/(nu-m-1)</span>
<a name="cl-77" style="color:rgb(53,114,176)" class=""></a> <span class="">RNGScope</span> <span class="">rngscope</span><span class="">;</span>
<a name="cl-78" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-79" style="color:rgb(53,114,176)" class=""></a> <span class="" style="color:rgb(68,85,136);font-weight:bold">int</span> <span class="">m</span> <span class="" style="font-weight:bold">=</span> <span class="">V</span><span class="">.</span><span class="">nrow</span><span class="">();</span>
<a name="cl-80" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-81" style="color:rgb(53,114,176)" class=""></a> <span class="">mat</span> <span class="">Vm</span><span class="">(</span><span class="">V</span><span class="">.</span><span class="">begin</span><span class="">(),</span> <span class="">m</span><span class="">,</span> <span class="">m</span><span class="">,</span> <span class="" style="color:rgb(153,153,153)">false</span><span class="">);</span>
<a name="cl-82" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-83" style="color:rgb(53,114,176)" class=""></a> <span class="" style="color:rgb(153,153,136);font-style:italic">// Can't vectorise because Rcpp-sugar rchisq doesnt vectorise df argument</span>
<a name="cl-84" style="color:rgb(53,114,176)" class=""></a> <span class="" style="color:rgb(153,153,136);font-style:italic">// T has sqrt chisqs on diagonal and normals below diagonal and garbage above diagonal</span>
<a name="cl-85" style="color:rgb(53,114,176)" class=""></a> <span class="">mat</span> <span class="">T</span><span class="">(</span><span class="">m</span><span class="">,</span><span class="">m</span><span class="">);</span>
<a name="cl-86" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-87" style="color:rgb(53,114,176)" class=""></a> <span class="" style="font-weight:bold">for</span><span class="">(</span><span class="" style="color:rgb(68,85,136);font-weight:bold">int</span> <span class="">i</span> <span class="" style="font-weight:bold">=</span> <span class="" style="color:rgb(0,153,153)">0</span><span class="">;</span> <span class="">i</span> <span class="" style="font-weight:bold"><</span> <span class="">m</span><span class="">;</span> <span class="">i</span><span class="" style="font-weight:bold">++</span><span class="">)</span> <span class="">{</span>
<a name="cl-88" style="color:rgb(53,114,176)" class=""></a> <span class="">T</span><span class="">(</span><span class="">i</span><span class="">,</span><span class="">i</span><span class="">)</span> <span class="" style="font-weight:bold">=</span> <span class="">sqrt</span><span class="">(</span><span class="">R</span><span class="" style="font-weight:bold">::</span><span class="">rchisq</span><span class="">(</span><span class="">nu</span><span class="" style="font-weight:bold">-</span><span class="">i</span><span class="">));</span>
<a name="cl-89" style="color:rgb(53,114,176)" class=""></a> <span class="">}</span>
<a name="cl-90" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-91" style="color:rgb(53,114,176)" class=""></a> <span class="" style="font-weight:bold">for</span><span class="">(</span><span class="" style="color:rgb(68,85,136);font-weight:bold">int</span> <span class="">j</span> <span class="" style="font-weight:bold">=</span> <span class="" style="color:rgb(0,153,153)">0</span><span class="">;</span> <span class="">j</span> <span class="" style="font-weight:bold"><</span> <span class="">m</span><span class="">;</span> <span class="">j</span><span class="" style="font-weight:bold">++</span><span class="">)</span> <span class="">{</span>
<a name="cl-92" style="color:rgb(53,114,176)" class=""></a> <span class="" style="font-weight:bold">for</span><span class="">(</span><span class="" style="color:rgb(68,85,136);font-weight:bold">int</span> <span class="">i</span> <span class="" style="font-weight:bold">=</span> <span class="">j</span><span class="" style="font-weight:bold">+</span><span class="" style="color:rgb(0,153,153)">1</span><span class="">;</span> <span class="">i</span> <span class="" style="font-weight:bold"><</span> <span class="">m</span><span class="">;</span> <span class="">i</span><span class="" style="font-weight:bold">++</span><span class="">)</span> <span class="">{</span>
<a name="cl-93" style="color:rgb(53,114,176)" class=""></a> <span class="">T</span><span class="">(</span><span class="">i</span><span class="">,</span><span class="">j</span><span class="">)</span> <span class="" style="font-weight:bold">=</span> <span class="">R</span><span class="" style="font-weight:bold">::</span><span class="">norm_rand</span><span class="">();</span>
<a name="cl-94" style="color:rgb(53,114,176)" class=""></a> <span class="">}}</span>
<a name="cl-95" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-96" style="color:rgb(53,114,176)" class=""></a> <span class="" style="color:rgb(153,153,136);font-style:italic">// explicitly declare T as triangular</span>
<a name="cl-97" style="color:rgb(53,114,176)" class=""></a> <span class="" style="color:rgb(153,153,136);font-style:italic">// top triangular is NaN ###</span>
<a name="cl-98" style="color:rgb(53,114,176)" class=""></a> <span class="">mat</span> <span class="">C</span> <span class="" style="font-weight:bold">=</span> <span class="">trans</span><span class="">(</span><span class="">trimatl</span><span class="">(</span><span class="">T</span><span class="">))</span> <span class="" style="font-weight:bold">*</span> <span class="">chol</span><span class="">(</span><span class="">Vm</span><span class="">);</span>
<a name="cl-99" style="color:rgb(53,114,176)" class=""></a> <span class="">mat</span> <span class="">CI</span> <span class="" style="font-weight:bold">=</span> <span class="">C</span><span class="">.</span><span class="">i</span><span class="">();</span>
<a name="cl-100" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-101" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// C is the upper triangular root of Wishart </span>
<a name="cl-102" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// therefore, W=C'C this is the LU decomposition </span>
<a name="cl-103" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// Inv(W) = CICI' this is the UL decomp *not LU*!</span>
<a name="cl-104" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-105" style="color:rgb(53,114,176)" class=""></a><span class="" style="color:rgb(153,153,136);font-style:italic">// W is Wishart draw, IW is W^-1</span>
<a name="cl-106" style="color:rgb(53,114,176)" class=""></a>
<a name="cl-107" style="color:rgb(53,114,176)" class=""></a> <span class="" style="font-weight:bold">return</span> <span class="">List</span><span class="" style="font-weight:bold">::</span><span class="">create</span><span class="">(</span>
<a name="cl-108" style="color:rgb(53,114,176)" class=""></a> <span class="">_</span><span class="">[</span><span class="" style="color:rgb(187,136,68)">"W"</span><span class="">]</span> <span class="" style="font-weight:bold">=</span> <span class="">trans</span><span class="">(</span><span class="">C</span><span class="">)</span> <span class="" style="font-weight:bold">*</span> <span class="">C</span><span class="">,</span>
<a name="cl-109" style="color:rgb(53,114,176)" class=""></a> <span class="">_</span><span class="">[</span><span class="" style="color:rgb(187,136,68)">"IW"</span><span class="">]</span> <span class="" style="font-weight:bold">=</span> <span class="">CI</span> <span class="" style="font-weight:bold">*</span> <span class="">trans</span><span class="">(</span><span class="">CI</span><span class="">),</span>
<a name="cl-110" style="color:rgb(53,114,176)" class=""></a> <span class="">_</span><span class="">[</span><span class="" style="color:rgb(187,136,68)">"C"</span><span class="">]</span> <span class="" style="font-weight:bold">=</span> <span class="">C</span><span class="">,</span>
<a name="cl-111" style="color:rgb(53,114,176)" class=""></a> <span class="">_</span><span class="">[</span><span class="" style="color:rgb(187,136,68)">"CI"</span><span class="">]</span> <span class="" style="font-weight:bold">=</span> <span class="">CI</span>
<a name="cl-112" style="color:rgb(53,114,176)" class=""></a> <span class="">);</span>
<a name="cl-113" style="color:rgb(53,114,176)" class=""></a><span class="">}</span></pre></div></div><div class="gmail_extra"><br class=""><div class="gmail_quote">On Thu, Jun 11, 2015 at 12:59 PM, Yue Li <span dir="ltr" class=""><<a href="mailto:gorillayue@gmail.com" target="_blank" class="">gorillayue@gmail.com</a>></span> wrote:<br class=""><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word" class="">Dear List,<div class=""><br class=""></div><div class="">I wonder if anyone could share their Rcpp code for an equivalence of rWishart. I’ve come across a similar question in the forum:</div><div class=""><br class=""></div><div class=""><a href="http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp" target="_blank" class="">http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp</a></div><div class=""><br class=""></div><div class="">Based on Dirk’s suggestion, I have my following attempt:</div><div class=""><br class=""></div><div class=""><div class="">// [[Rcpp::export]]</div><div class="">mat mvrnormArma(int n, vec mu, mat sigma) {</div><div class=""> </div><div class=""> int ncols = sigma.n_cols;</div><div class=""> </div><div class=""> mat Y = randn(n, ncols);</div><div class=""> </div><div class=""> return repmat(mu, 1, n).t() + Y * chol(sigma);</div><div class="">}</div><div class=""><br class=""></div><div class=""><br class=""></div><div class="">// [[Rcpp::export]]</div><div class="">mat rWishartCpp(mat sigma, int n) {</div><div class=""> </div><div class=""> mat X = mvrnormArma(1, zeros<vec>(sigma.n_rows), sigma);</div><div class=""> </div><div class=""> return X.t() * X;</div><div class="">}</div></div><div class=""><br class=""></div><div class="">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:</div><div class=""><br class=""></div><div class="">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 …</div><div class=""><br class=""></div><div class=""><div class="">function baysm_rwishart(nu, V)</div><div class="">{</div><div class=""> m = nrow(V)</div><div class=""> </div><div class=""> df = (nu + nu - m + 1) - (nu - m + 1):nu</div><div class=""> </div><div class=""> if (m > 1) {</div><div class=""> T = diag(sqrt(rchisq(c(rep(1, m)), df)))</div><div class=""> T[lower.tri(T)] = rnorm((m * (m + 1)/2 - m))</div><div class=""> } </div><div class=""> else {</div><div class=""> T = sqrt(rchisq(1, df))</div><div class=""> }</div><div class=""> </div><div class=""> U = chol(V)</div><div class=""> </div><div class=""> C = t(T) %*% U</div><div class=""> </div><div class=""> CI = backsolve(C, diag(m))</div><div class=""> </div><div class=""> return(list(W = crossprod(C), </div><div class=""> <span style="white-space:pre-wrap" class=""> </span>IW = crossprod(t(CI)), C = C, </div><div class=""> CI = CI))</div><div class="">}</div></div><span class="HOEnZb"><font color="#888888" class=""><div class=""><br class=""></div><div class="">Yue</div><div class=""><br class=""></div><div class=""><br class=""></div></font></span></div><br class="">_______________________________________________<br class="">
Rcpp-devel mailing list<br class="">
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" class="">Rcpp-devel@lists.r-forge.r-project.org</a><br class="">
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank" class="">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br class=""></blockquote></div><br class=""></div>
</div></blockquote></div><br class=""></div></body></html>