[R] Need help to understand integrate function

[R_help Help]

Hi Ravi - Thank you. I'm wondering how the magic happens here with
options(digits=10)? And the most important point, I do not quite
understand why divergent error could occur. In this case, the function
is analytically integrable. Any insight would be greatly appreciated
so that I learn something out of this. Thank you.

Rob

On Wed, Aug 11, 2010 at 11:46 PM, Ravi Varadhan <rvaradhan at jhmi.edu> wrote:
> Try this:
>
> options(digits=10)
>
>> integrate(f=powerLaw2, lower=0, upper=Inf, l1=1.8980185, l2=-0.0804259, c0=1, t0=259.78, rel.tol=1.e-10)
> 0.01089019946 with absolute error < 3.7e-11
>
>
> Ravi.
> ____________________________________________________________________
>
> Ravi Varadhan, Ph.D.
> Assistant Professor,
> Division of Geriatric Medicine and Gerontology
> School of Medicine
> Johns Hopkins University
>
> Ph. (410) 502-2619
> email: rvaradhan at jhmi.edu
>
>
> ----- Original Message -----
> From: R_help Help <rhelpacc at gmail.com>
> Date: Wednesday, August 11, 2010 9:44 pm
> Subject: [R] Need help to understand integrate function
> To: r-help at r-project.org
>
>
>> Hi,
>>
>>  I'm running into a wall when trying to use the integrate function. I
>>  have the following setting:
>>
>>  powerLaw2 <- function(x,l1,l2,c0,t0) {
>>
>>     idx <- which(x <= 0);
>>     if (length(idx) > 0) {
>>        x[idx] <- 0;
>>     }
>>
>>     xl <- (-l1+l2)*log(x/t0);
>>     L <- log(c0)-l1*log(x)-log(1+exp(xl));
>>     L <- exp(L);
>>     return(L);
>>  }
>>
>>  plCDF2 <- function(x,l1,l2,c0,t0) {
>>
>>     print(c(l1,l2,c0,t0));
>>     cdf <- integrate(f=powerLaw2,lower=x,upper=Inf,l1=l1,l2=l2,c0=c0,t0=t0,stop.on.error=T);
>>     return(cdf$value);
>>
>>  }
>>
>>  I know that as long as my l1 > 1 and l2 < 1 the integration above
>>  should converge. However, this doesn't seem to be the case when I call
>>  integrate function with the following parameter:
>>
>>
>>  Browse[1]> integrate(f=powerLaw2,lower=1e-09,upper=Inf,l1=1.8980185,l2=-0.0804259,c0=1,t0=259.78,stop.on.error=T)
>>  failed with message 'the integral is probably divergent'
>>
>>  Would you please shed some light? And how could I prevent this
>>  situation? Thank you.
>>
>>  Rgds,
>>
>>  Rob
>>
>>  ______________________________________________
>>  R-help at r-project.org mailing list
>>
>>  PLEASE do read the posting guide
>>  and provide commented, minimal, self-contained, reproducible code.
>