<div dir="ltr"><div><div><div>Dear Dirk, thanks for your answer, I'll keep studying and decide what way to follow.<br><br>Dear Qiang, here's the code for the function (I already did some work to the math behind the density itself in order to quicken a bit)<br><br>library(cubature)<br>library(zipfR)<br><br>dcbinom <- Vectorize(function(x,nu,pp) { # density function<br> <br> stopifnot(pp>0&pp<1&x<=nu+1&nu>0)<br> if (x<=0|x>=nu+1) 0 <br> else {<br> ifelse(x>(nu+1)/2&x<nu&pp>.5, {<br> aa <- exp(-2*lbeta(x,nu-x+1))<br> bb <- function(s) (s[1]*s[2])^(nu-x)*((1-s[1])*(1-s[2]))^(x-1)*(log((1-s[2])/s[2])-log((1-s[1])/s[1])) <br> aa*adaptIntegrate(bb,c(0,0),c(1,1-pp))$integral<br> },<br>{<br> aa <- exp(-2*lbeta(x,nu-x+1))<br> bb <- function(s) (s[1]*s[2])^(x-1)*((1-s[1])*(1-s[2]))^(nu-x)*(log((1-s[2])/s[2])-log((1-s[1])/s[1])) <br> aa*adaptIntegrate(bb,c(pp,1E-10),c(1-1E-10,pp))$integral<br> <br>})}<br>})<br><br><br></div>Now, this is *quite* fast if I want to obtain a single density value. Still, if I need to use it repeatedly to, e.g., draw the density using curve(), it gets quite slow:<br><br><span class="" style="border-collapse:separate;color:rgb(0,0,0);font-family:"Lucida Console";font-size:13px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:15px;text-indent:0px;text-transform:none;white-space:pre-wrap;word-spacing:0px;background-color:rgb(225,226,229)"><pre tabindex="0" class="" id="rstudio_console_output" style="font-family:"Lucida Console";font-size:10pt!important;outline-style:none;outline-width:initial;outline-color:initial;border-style:none;border-width:initial;border-color:initial;white-space:pre-wrap!important;word-break:break-all;margin:0px;line-height:1.2"><span class="" style="color:blue">> system.time(curve(dcbinom(x,nu=10,pp=.2),xlim=c(0,11)))
</span> user system elapsed
1.48 0.00 1.48 </pre></span><br><span class="" style="border-collapse:separate;color:rgb(0,0,0);font-family:"Lucida Console";font-size:13px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:15px;text-indent:0px;text-transform:none;white-space:pre-wrap;word-spacing:0px;background-color:rgb(225,226,229)"><pre tabindex="0" class="" id="rstudio_console_output" style="font-family:"Lucida Console";font-size:10pt!important;outline-style:none;outline-width:initial;outline-color:initial;border-style:none;border-width:initial;border-color:initial;white-space:pre-wrap!important;word-break:break-all;margin:0px;line-height:1.2"><span class="" style="white-space:pre;color:blue">> </span><span class="" style="color:blue">system.time(curve(dcbinom(x,nu=1,pp=.7),xlim=c(0,2)))
</span> user system elapsed
4.25 0.00 4.25 <br><br></pre><pre tabindex="0" class="" id="rstudio_console_output" style="font-family:"Lucida Console";font-size:10pt!important;outline-style:none;outline-width:initial;outline-color:initial;border-style:none;border-width:initial;border-color:initial;white-space:pre-wrap!important;word-break:break-all;margin:0px;line-height:1.2">Moreover, I have the feeling that something's wrong with the numerical values when approaching the limits of the support, i.e. (0, nu+1). <br></pre></span><br></div>Any suggestion?<br><br></div>Thanks again!<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jan 9, 2015 at 3:42 PM, Qiang Kou <span dir="ltr"><<a href="mailto:qkou@umail.iu.edu" target="_blank">qkou@umail.iu.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi,<div><br></div><div>The backend of cubature package is the cubature C library from MIT. The same author with nlopt.</div><div><br></div><div>Can you give a little more detail about the double integral? Cubature should not be that slow from my experience.</div><div><br></div><div>Best,</div><div><br></div><div>KK</div><div><br></div></div><div class="gmail_extra"><div><div class="h5"><br><div class="gmail_quote">On Fri, Jan 9, 2015 at 9:38 AM, Dirk Eddelbuettel <span dir="ltr"><<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div><br>
On 9 January 2015 at 15:10, Federico Andreis wrote:<br>
| Dear all,<br>
|<br>
| I was wondering what, in your opinion, would be the best way to handle the<br>
| computation of definite integrals using Rcpp.<br>
|<br>
| To be more precise, I need to evaluate a double integral that turns out to be<br>
| somewhat nasty (really slow computation and presumably inaccurate results using<br>
| the R cubature package).<br>
|<br>
| Should I write the integration algorithm from scratch, or is there any external<br>
| library you would suggest?<br>
|<br>
| I've already found the post on Stackoverflow 'using C function from other<br>
| package in Rcpp' that dealt with an integration problem as well, but it looked<br>
| too general for my problem and also, it's one year old, maybe something else<br>
| has turned up in the meanwhile..<br>
<br>
</div></div>From the top of my head, I think there are several CRAN package doing<br>
something like cubature with related names. It is not something I use day<br>
to day so maybe somebody else will chime in...<br>
<br>
In general, though, one (easy ?) approach is to solve the problem with<br>
standalone C/C++ code [ possibly entirely outside of of R ] and then learn<br>
(here and elsewhere) how to connect that to R via Rcpp. It is<br>
straightforward.<br>
<br>
Alternatively, and starting from the other end, learn how to access C / C++<br>
code in another package and do it with C++ code from R via Rcpp.<br>
<br>
Or use the middle ground. There are lots of ways in which Rcpp can help you here.<br>
<br>
Dirk<br>
<span><br>
| Thanks in advance, and congrats for the great work with Rcpp!<br>
|<br>
| /federico<br>
</span>| _______________________________________________<br>
| Rcpp-devel mailing list<br>
| <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">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>
<span><font color="#888888"><br>
--<br>
<a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a> | @eddelbuettel | <a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a><br>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">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>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br></div></div><span class="HOEnZb"><font color="#888888"><div><div dir="ltr">Qiang Kou<div><a href="mailto:qkou@umail.iu.edu" target="_blank">qkou@umail.iu.edu</a><br><div>School of Informatics and Computing, Indiana University</div><div><br></div></div></div></div>
</font></span></div>
</blockquote></div><br></div>