[R] integrate with large parameters

[Duncan Murdoch]

On 01/05/2009 5:46 AM, Andreas Wittmann wrote:
> Dear R-users,
> 
> i have to integrate the following function
> 
> `fun1` <- function (a, l1, l2)
> {
>   exp(log(l1) * (a - 1) - l2 * lgamma(a))
> }
> 
> but if l1 is large, i get the "non-finite function value" error, so my 
> idea is to rescale with exp(-l1)
> 
> `fun2` <- function (a, l1, l2)
> {
>   exp(log(l1) * (a - 1) - l2 * lgamma(a) - l1)
> }
> 
> but it seems this doesn't solve the problem, when l1=1.0e80. maybe i 
> have to iterate between large values in integrate and large values in 
> exp(l1) and then use the normal or the rescaled version?

I think if you can work out the mode of the function, and divide by 
that, it should be possible to make things work.  (The main problem will 
be that a large range of function values will underflow to zero. 
Hopefully those have a negligible integral, but if integrate evaluates 
the target there too often, it might think the integral is zero.)

Duncan Murdoch

> 
> 
>  > l1<-10; l2<-10;
>  > integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
> [1] 11.65311
>  > integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val 
> * exp(l1)
> [1] 11.65312
>  >
>  > l1<-1.0e12; l2=50;
>  > integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
> [1] 929558588152
>  > integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val 
> * exp(l1)
> [1] NaN
>  >
>  > l1<-1.0e20; l2=50;
>  > integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
> [1] 5.005192e+24
>  > integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val 
> * exp(l1)
> [1] NaN
>  >
>  > l1<-1.0e80; l2=50;
>  > integrate(Vectorize(function(a) {fun1(a, l1=l1, l2=l2)}), 1, Inf)$val
> Fehler in integrate(Vectorize(function(a) { : non-finite function value
>  > integrate(Vectorize(function(a) {fun2(a, l1=l1, l2=l2)}), 1, Inf)$val 
> * exp(l1)
> [1] NaN
> 
> 
> Many thanks if you have any advice for me!
> 
> best regards
> 
> Andreas
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.