[R] Improving result from integrate

[. .]

So guys,

I am writing to ask for an opinion from you...

As I begin to write in the previous post, I developed a function
called myintegrate that is based in integrate but my function give
better results for the problem I am dealing. My development was
inspired in this post
http://r.789695.n4.nabble.com/integration-td843252.html, by Ravi
Varadhan.

To check if the results are good I compare with the ones from an
analytical version of the problem.

This is the new code,

dsad <- Vectorize(FUN=
 function(y, frac, rate, sad, samp="Poisson", trunc=0, ...){
   f0 <- function(y,frac,n) {
     f1 <- function(y,frac,n){
       dpois(y,frac*n)
     }
     dcom <- paste("d",deparse(substitute(sad)),sep="")
     dots <- c(as.name("n"),list(...))
     f2 <- call(dcom,dots)
     f <- function(n){
       f1(y,frac,n)*f2(n,rate)
     }
     myintegrate <- function() {
       r <- 0
       r1 <- 1
       x1 <- 0
       dx <- 20
       while(r1 > 10e-500){
         r1 <- integrate(f,x1,x1+dx)$value
         r <- r + r1
         x1 <- x1 + dx
       }
       r + integrate(f,x1,Inf)$valu
     }
     myintegrate()
   }
   f0(y,frac,n)/(1-f0(trunc,frac,n))
 },"y")
dsad(200, 0.1, 0.1, exp)

whose result is 6.223015e-61, correct.

And the below is the old one,

dsad <- Vectorize(FUN=
 function(y, frac, rate, sad, samp="Poisson", trunc=0, ...){
   f0 <- function(y,frac,n) {
     f1 <- function(y,frac,n){
       dpois(y,frac*n)
     }
     dcom <- paste("d",deparse(substitute(sad)),sep="")
     dots <- c(as.name("n"),list(...))
     f2 <- call(dcom,dots)
     f <- function(n){
       f1(y,frac,n)*f2(n,rate)
     }
     myintegrate <- function() {
       r <- 0
       r1 <- 1
       x1 <- 0
       integrate(f,x1,Inf)$valu
     }
     myintegrate()
   }
   f0(y,frac,n)/(1-f0(trunc,frac,n))
 },"y")
dsad(200, 0.1, 0.1, exp)

whose result is 6.554036e-80, incorrect.

So, what do you think of it?

I will really apprettiate any comments.

Thanks in advance.

On Mon, Aug 29, 2011 at 4:06 AM, . . <xkziloj at gmail.com> wrote:
> Hi all,
>
> I am utilizing integrate() to get the area below a function that shape
> like an exponential, getting very low for higher values of the
> horizontal coordinate y.
>
> In the code below, I present the tricky way I can improve the result
> obtained. I compare this result with the one I get from Mathematica
> and also with the exact solution I have for this particular case.
>
> /*----------
> func <- function(y, a, rate){
>  x <- function(n,rate) {
>    rate*exp(-n*rate)
>  }
>  boi <- function(y,n,a){
>                      w <- y*log(a*n)-lfactorial(y)-a*n
>                      exp(w)
>                    }
>  f <- function(n){
>    boi(y,n,a)*x(n,rate)
>  }
>  r <- 0
>                    r1 <- 1
>                    x1 <- 0
>                    dx <- 20
>                    while(r1 > 10e-1000){
>                      r1 <- integrate(f,x1,x1+dx)$value
>                      r <- r + r1
>                      x1 <- x1 + dx
>                    }
>                    r + integrate(f,x1,Inf)$valu
> }
> func(200,0.1,0.1)
> ----------*/
>
> Altought I get better results, the value of dx must be carefully
> selected. So I ask, there is another method that can give me better
> results?
>