[Rcpp-devel] Integrals with Rcpp
Federico Andreis
federico.andreis at gmail.com
Fri Jan 9 16:26:36 CET 2015
Dear Dirk, thanks for your answer, I'll keep studying and decide what way
to follow.
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)
library(cubature)
library(zipfR)
dcbinom <- Vectorize(function(x,nu,pp) { # density function
stopifnot(pp>0&pp<1&x<=nu+1&nu>0)
if (x<=0|x>=nu+1) 0
else {
ifelse(x>(nu+1)/2&x<nu&pp>.5, {
aa <- exp(-2*lbeta(x,nu-x+1))
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]))
aa*adaptIntegrate(bb,c(0,0),c(1,1-pp))$integral
},
{
aa <- exp(-2*lbeta(x,nu-x+1))
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]))
aa*adaptIntegrate(bb,c(pp,1E-10),c(1-1E-10,pp))$integral
})}
})
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:
> system.time(curve(dcbinom(x,nu=10,pp=.2),xlim=c(0,11))) user system elapsed
1.48 0.00 1.48
> system.time(curve(dcbinom(x,nu=1,pp=.7),xlim=c(0,2))) user system elapsed
4.25 0.00 4.25
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).
Any suggestion?
Thanks again!
On Fri, Jan 9, 2015 at 3:42 PM, Qiang Kou <qkou at umail.iu.edu> wrote:
> Hi,
>
> The backend of cubature package is the cubature C library from MIT. The
> same author with nlopt.
>
> Can you give a little more detail about the double integral? Cubature
> should not be that slow from my experience.
>
> Best,
>
> KK
>
>
> On Fri, Jan 9, 2015 at 9:38 AM, Dirk Eddelbuettel <edd at debian.org> wrote:
>
>>
>> On 9 January 2015 at 15:10, Federico Andreis wrote:
>> | Dear all,
>> |
>> | I was wondering what, in your opinion, would be the best way to handle
>> the
>> | computation of definite integrals using Rcpp.
>> |
>> | To be more precise, I need to evaluate a double integral that turns out
>> to be
>> | somewhat nasty (really slow computation and presumably inaccurate
>> results using
>> | the R cubature package).
>> |
>> | Should I write the integration algorithm from scratch, or is there any
>> external
>> | library you would suggest?
>> |
>> | I've already found the post on Stackoverflow 'using C function from
>> other
>> | package in Rcpp' that dealt with an integration problem as well, but it
>> looked
>> | too general for my problem and also, it's one year old, maybe something
>> else
>> | has turned up in the meanwhile..
>>
>> From the top of my head, I think there are several CRAN package doing
>> something like cubature with related names. It is not something I use day
>> to day so maybe somebody else will chime in...
>>
>> In general, though, one (easy ?) approach is to solve the problem with
>> standalone C/C++ code [ possibly entirely outside of of R ] and then learn
>> (here and elsewhere) how to connect that to R via Rcpp. It is
>> straightforward.
>>
>> Alternatively, and starting from the other end, learn how to access C /
>> C++
>> code in another package and do it with C++ code from R via Rcpp.
>>
>> Or use the middle ground. There are lots of ways in which Rcpp can help
>> you here.
>>
>> Dirk
>>
>> | Thanks in advance, and congrats for the great work with Rcpp!
>> |
>> | /federico
>> | _______________________________________________
>> | Rcpp-devel mailing list
>> | Rcpp-devel at lists.r-forge.r-project.org
>> | https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>>
>> --
>> http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org
>> _______________________________________________
>> Rcpp-devel mailing list
>> Rcpp-devel at lists.r-forge.r-project.org
>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>>
>
>
>
> --
> Qiang Kou
> qkou at umail.iu.edu
> School of Informatics and Computing, Indiana University
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20150109/0ce285d5/attachment-0001.html>
More information about the Rcpp-devel
mailing list