<html dir="ltr">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=Windows-1252">
<style>
<!--
@font-face
{font-family:Calibri}
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif"}
a:link, span.MsoHyperlink
{color:blue;
text-decoration:underline}
a:visited, span.MsoHyperlinkFollowed
{color:purple;
text-decoration:underline}
span.Shkpostityyli17
{font-family:"Calibri","sans-serif";
color:windowtext}
span.keyword
{}
span.keywordtype
{}
span.lineno
{}
span.SpellE
{}
.MsoChpDefault
{font-family:"Calibri","sans-serif"}
@page WordSection1
{margin:72.0pt 72.0pt 72.0pt 72.0pt}
div.WordSection1
{}
-->
</style><style id="owaParaStyle" type="text/css">P {margin-top:0;margin-bottom:0;}</style>
</head>
<body ocsi="0" fpstyle="1" style="" lang="EN-US" link="blue" vlink="purple">
<div style="direction: ltr;font-family: Tahoma;color: #000000;font-size: 10pt;">
<div class="WordSection1">
<p class="MsoNormal">Dear list,</p>
<p class="MsoNormal"> </p>
<p class="MsoNormal">I am trying to call R’s negative binomial function with <span class="SpellE">
Rcpp</span>, but encountered some weird behaviour. Here is an example using inline:</p>
<p class="MsoNormal"> </p>
<br>
rcpp_Rf_dnbinom <- rcpp(signature(),<br>
' return wrap(Rf_dnbinom( 4.0, 0.5, 0.9, 1)); ')<br>
<br>
rcpp_Rf_dnbinom_mu <- rcpp(signature(),<br>
' return wrap(Rf_dnbinom_mu( 4.0, 0.5, 0.9, 1)); ')<br>
<br>
<br>
rcpp_dnbinom <- rcpp(signature(),<br>
' return wrap(R::dnbinom( 4.0, 0.5, 0.9, 1)); ')<br>
<br>
rcpp_dnbinom_mu <- rcpp(signature(),<br>
' return wrap(R::dnbinom_mu( 4.0, 0.5, 0.9, 1)); ')<br>
<br>
rcpp_dnbinom_sugar <- rcpp(signature(y="numeric"), <br>
' NumericVector x = NumericVector(y);<br>
NumericVector res = dnbinom( x, 0.5, 0.9,1);<br>
return wrap(res); ')<br>
<br>
<br>
rcpp_dnbinom_mu_sugar <- rcpp(signature(y="numeric"), <br>
' NumericVector x = NumericVector(y);<br>
NumericVector res = dnbinom_mu( x, 0.5, 0.9,1);<br>
return wrap(res); ')<br>
<br>
rcpp_Rf_dnbinom()<br>
#-10.5597<br>
rcpp_Rf_dnbinom_mu()<br>
#-3.578823<br>
rcpp_dnbinom()<br>
#-10.5597<br>
rcpp_dnbinom_mu()<br>
#-10.5597<br>
rcpp_dnbinom_sugar(y=4.0)<br>
#-10.5597<br>
rcpp_dnbinom_mu_sugar(y=4.0)<br>
#-3.578823<br>
dnbinom(x=4,size=0.5,mu=0.9,log=TRUE)<br>
#[1] -3.578823<br>
dnbinom(x=4,size=0.5,prob=0.9,log=TRUE)<br>
#[1] -10.5597<br>
<p class="MsoNormal"> </p>
<p class="MsoNormal">So it looks like that everything is fine when using <span class="SpellE">
Rcpp</span> sugar or Rf_dbinom_mu directly, but when using form R::dbinom both <span class="SpellE">
dnbinom</span> and <span class="SpellE">dnbinom_mu</span> <span style=""> </span>calls actually use
<span class="SpellE">dnbinom</span>. </p>
<p class="MsoNormal"><br>
</p>
<p class="MsoNormal"><span style=""> </span>I am not sure if this is relevant, but
<span class="SpellE">Rmath.h</span> in Rcpp contains lines:</p>
<p class="MsoNormal"> </p>
<p class="MsoNormal"><span class="keyword">129 inline</span> <span class="keywordtype">
double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a1bdf703fb4850bb68382cef265fdc0c6" target="_blank">
<span class="SpellE">dnbinom</span></a>(<span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10" target="_blank">
x</a>, <span class="keywordtype">double</span> <span class="SpellE">sz</span>, <span class="keywordtype">
double</span> <span class="SpellE">pb</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lg</span>) { return ::<span class="SpellE">Rf_dnbinom</span>(x,
<span class="SpellE">sz</span>, <span class="SpellE">pb</span>, <span class="SpellE">
lg</span>); }</p>
<p class="MsoNormal"><a name="l00130"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a2e59e21ed1007c3cc9394b39d0c04e1d" target="_blank">130</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a2e59e21ed1007c3cc9394b39d0c04e1d" target="_blank">
<span class="SpellE">pnbinom</span></a>(<span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10" target="_blank">
x</a>, <span class="keywordtype">double</span> <span class="SpellE">sz</span>, <span class="keywordtype">
double</span> <span class="SpellE">pb</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lt</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lg</span>) { return ::<span class="SpellE">Rf_pnbinom</span>(x,
<span class="SpellE">sz</span>, <span class="SpellE">pb</span>, <span class="SpellE">
lt</span>, <span class="SpellE">lg</span>); }</p>
<p class="MsoNormal"><a name="l00131"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a6152e2d7c265f629acdf0adae6a90989" target="_blank">131</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a6152e2d7c265f629acdf0adae6a90989" target="_blank">
<span class="SpellE">qnbinom</span></a>(<span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/external__pointer_8r.html#a745dfbf3bbf4ccff97d7b764f8694d25" target="_blank">
p</a>, <span class="keywordtype">double</span> <span class="SpellE">sz</span>, <span class="keywordtype">
double</span> <span class="SpellE">pb</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lt</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lg</span>) { return ::<span class="SpellE">Rf_qnbinom</span>(p,
<span class="SpellE">sz</span>, <span class="SpellE">pb</span>, <span class="SpellE">
lt</span>, <span class="SpellE">lg</span>); }</p>
<p class="MsoNormal"><a name="l00132"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ad955db033a0d939ed1b37d454a3a6b29" target="_blank">132</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ad955db033a0d939ed1b37d454a3a6b29" target="_blank">
<span class="SpellE">rnbinom</span></a>(<span class="keywordtype">double</span> <span class="SpellE">
sz</span>, <span class="keywordtype">double</span> <span class="SpellE">pb</span>) { return ::<span class="SpellE">Rf_rnbinom</span>(<span class="SpellE">sz</span>,
<span class="SpellE">pb</span>); }</p>
<p class="MsoNormal"><a name="l00133"></a><span class="lineno">133</span> </p>
<p class="MsoNormal"><a name="l00134"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a351c1ce1012d1a965edaaee71cfb4031" target="_blank">134</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a351c1ce1012d1a965edaaee71cfb4031" target="_blank">
<span class="SpellE">dnbinom_mu</span></a>(<span class="keywordtype">double</span>
<a href="http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10" target="_blank">
x</a>, <span class="keywordtype">double</span> <span class="SpellE">sz</span>, <span class="keywordtype">
double</span> mu, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lg</span>) { return ::<span class="SpellE">Rf_dnbinom</span>(x,
<span class="SpellE">sz</span>, mu, <span class="SpellE">lg</span>); }</p>
<p class="MsoNormal"><a name="l00135"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#aea24f458a776074cba6a9e2ecffaecc6" target="_blank">135</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#aea24f458a776074cba6a9e2ecffaecc6" target="_blank">
<span class="SpellE">pnbinom_mu</span></a>(<span class="keywordtype">double</span>
<a href="http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10" target="_blank">
x</a>, <span class="keywordtype">double</span> <span class="SpellE">sz</span>, <span class="keywordtype">
double</span> mu, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lt</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lg</span>) { return ::<span class="SpellE">Rf_pnbinom</span>(x,
<span class="SpellE">sz</span>, mu, <span class="SpellE">lt</span>, <span class="SpellE">
lg</span>); }</p>
<p class="MsoNormal"><a name="l00136"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a5e4cb4981198b228e9791afb93caed4e" target="_blank">136</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a5e4cb4981198b228e9791afb93caed4e" target="_blank">
<span class="SpellE">qnbinom_mu</span></a>(<span class="keywordtype">double</span>
<a href="http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10" target="_blank">
x</a>, <span class="keywordtype">double</span> <span class="SpellE">sz</span>, <span class="keywordtype">
double</span> mu, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lt</span>, <span class="SpellE"><span class="keywordtype">int</span></span>
<span class="SpellE">lg</span>) { return ::<span class="SpellE">Rf_qnbinom</span>(x,
<span class="SpellE">sz</span>, mu, <span class="SpellE">lt</span>, <span class="SpellE">
lg</span>); }</p>
<p class="MsoNormal"><a name="l00137"></a><span class="lineno"><a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ae51cfffdc309fec3380b0e88c01ddc2e" target="_blank">137</a></span>
<span class="keyword">inline</span> <span class="keywordtype">double</span> <a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ae51cfffdc309fec3380b0e88c01ddc2e" target="_blank">
<span class="SpellE">rnbinom_mu</span></a>(<span class="keywordtype">double</span>
<span class="SpellE">sz</span>, <span class="keywordtype">double</span> mu) { return ::<span class="SpellE">Rf_rnbinom</span>(<span class="SpellE">sz</span>, mu); }</p>
<p class="MsoNormal"> </p>
<p class="MsoNormal">To me this <span class="SpellE">looks</span> like that both <span class="SpellE">
dnbinom</span> and <span class="SpellE">dnbinom_mu</span> calls the same function.</p>
<p class="MsoNormal"> </p>
<p class="MsoNormal">So, is there something wrong with my code or is there a typo in R::<a href="http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a351c1ce1012d1a965edaaee71cfb4031" target="_blank">
<span class="SpellE">dnbinom_mu</span></a>?</p>
<p class="MsoNormal"><br>
</p>
<p class="MsoNormal">To be honest, I don't really understand the differences between the three versions, except that the sugar version needs NumericVector as first argument and it is vectorized regards that parameter. I would actually need version which is
vectorized wrt all arguments (except log), but since there isn't one I am doing something like this:<br>
</p>
<p class="MsoNormal"><br>
</p>
<p class="MsoNormal">#include "RcppArmadillo.h"<br>
</p>
<p class="MsoNormal">// [[Rcpp::depends(RcppArmadillo)]]<br>
</p>
<p class="MsoNormal">double fun(const arma::mat& y, const arma::mat& x ,const arma::mat& theta){<br>
double res=0.0; <br>
for(unsigned int i=0; i<y.n_elem; i++){<br>
if(arma::is_finite(y(i))){<br>
res += Rf_dbinom_mu( y(i), u(i), theta(i), 1); //was R::dbinom_mu<br>
}<br>
}<br>
</p>
<p class="MsoNormal"> return res;</p>
<p class="MsoNormal">}</p>
<p class="MsoNormal"><br>
</p>
<p class="MsoNormal"><br>
</p>
<p class="MsoNormal">Best regards,</p>
<p class="MsoNormal"> </p>
<p class="MsoNormal">Jouni Helske</p>
<p class="MsoNormal"> </p>
</div>
</div>
</body>
</html>