<div dir="ltr">I had ported the bayesm version a while a back, here it is:<div><br></div><div><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)"><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)"></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)"></a><span class="" style="color:rgb(153,153,136);font-style:italic">// </span>
<a name="cl-40" style="color:rgb(53,114,176)"></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)"></a><span class="" style="color:rgb(153,153,136);font-style:italic">// </span>
<a name="cl-42" style="color:rgb(53,114,176)"></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)"></a> <span class="">RNGScope</span> <span class="">rngscope</span><span class="">;</span>
<a name="cl-78" style="color:rgb(53,114,176)"></a>
<a name="cl-79" style="color:rgb(53,114,176)"></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)"></a>
<a name="cl-81" style="color:rgb(53,114,176)"></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)"></a>
<a name="cl-83" style="color:rgb(53,114,176)"></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)"></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)"></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)"></a>
<a name="cl-87" style="color:rgb(53,114,176)"></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)"></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)"></a> <span class="">}</span>
<a name="cl-90" style="color:rgb(53,114,176)"></a>
<a name="cl-91" style="color:rgb(53,114,176)"></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)"></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)"></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)"></a> <span class="">}}</span>
<a name="cl-95" style="color:rgb(53,114,176)"></a>
<a name="cl-96" style="color:rgb(53,114,176)"></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)"></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)"></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)"></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)"></a>
<a name="cl-101" style="color:rgb(53,114,176)"></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)"></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)"></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)"></a>
<a name="cl-105" style="color:rgb(53,114,176)"></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)"></a>
<a name="cl-107" style="color:rgb(53,114,176)"></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)"></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)"></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)"></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)"></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)"></a> <span class="">);</span>
<a name="cl-113" style="color:rgb(53,114,176)"></a><span class="">}</span></pre></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Jun 11, 2015 at 12:59 PM, Yue Li <span dir="ltr"><<a href="mailto:gorillayue@gmail.com" target="_blank">gorillayue@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">Dear List,<div><br></div><div>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><br></div><div><a href="http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp" target="_blank">http://stackoverflow.com/questions/23463852/generate-wishart-distribution-using-rwishart-within-rcpp</a></div><div><br></div><div>Based on Dirk’s suggestion, I have my following attempt:</div><div><br></div><div><div>// [[Rcpp::export]]</div><div>mat mvrnormArma(int n, vec mu, mat sigma) {</div><div> </div><div> int ncols = sigma.n_cols;</div><div> </div><div> mat Y = randn(n, ncols);</div><div> </div><div> return repmat(mu, 1, n).t() + Y * chol(sigma);</div><div>}</div><div><br></div><div><br></div><div>// [[Rcpp::export]]</div><div>mat rWishartCpp(mat sigma, int n) {</div><div> </div><div> mat X = mvrnormArma(1, zeros<vec>(sigma.n_rows), sigma);</div><div> </div><div> return X.t() * X;</div><div>}</div></div><div><br></div><div>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><br></div><div>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><br></div><div><div>function baysm_rwishart(nu, V)</div><div>{</div><div> m = nrow(V)</div><div> </div><div> df = (nu + nu - m + 1) - (nu - m + 1):nu</div><div> </div><div> if (m > 1) {</div><div> T = diag(sqrt(rchisq(c(rep(1, m)), df)))</div><div> T[lower.tri(T)] = rnorm((m * (m + 1)/2 - m))</div><div> } </div><div> else {</div><div> T = sqrt(rchisq(1, df))</div><div> }</div><div> </div><div> U = chol(V)</div><div> </div><div> C = t(T) %*% U</div><div> </div><div> CI = backsolve(C, diag(m))</div><div> </div><div> return(list(W = crossprod(C), </div><div> <span style="white-space:pre-wrap"> </span>IW = crossprod(t(CI)), C = C, </div><div> CI = CI))</div><div>}</div></div><span class="HOEnZb"><font color="#888888"><div><br></div><div>Yue</div><div><br></div><div><br></div></font></span></div><br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></blockquote></div><br></div>