"non-finite function value" when using integrate() in R

Question

I am using the following R code, taken from a published paper (citation below). This is the code:

int2=function(x,r,n,p) {
    (1+x)^((n-1-p)/2)*(1+(1-r^2)*x)^(-(n-1)/2)*x^(-3/2)*exp(-n/(2*x))}
integrate(f=int2,lower=0,upper=Inf,n=530,r=sqrt(.245),p=3, stop.on.error=FALSE)

When I run it, I get the error "non-finite function value". Yet Maple is able to compute this as 4.046018765*10^27.

I tried using "integral" in package pracma, which gives me a different error: Error in if (delta < tol) break : missing value where TRUE/FALSE needed

The overall goal is to compute a ratio of two integrals, as described in Wetzels & Wagenmakers (2012) "A default Bayesian hypothesis test for correlations" (http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3505519/). The entire function is as follows:

jzs.pcorbf = function(r0, r1, p0, p1, n) {
  int = function(r,n,p,g) {
    (1+g)^((n-1-p)/2)*(1+(1-r^2)*g)^(-(n-1)/2)*g^(-3/2)*exp(-n/(2*g))};
  bf10=integrate(int, lower=0,upper=Inf,r=r1,p=p1,n=n)$value/
       integrate(int,lower=0,upper=Inf,r=r0,p=p0,n=n)$value;
  return(bf10)
}

Thanks!

I don't think that is it. There are no missing or NaN input values. There is no data. It's the integration of a simple mathematical expression. — Dimitri
That's not a good way of looking at it, Dimitri: numerical integration does use "data," as it has to calculate the value of the integrand at a specified number of locations. — Carl Witthoft

josliber josliber · Accepted Answer · 2015-02-22T03:06:01

The issue is that your integral function is generating NaN values when called with x values in its domain. You're integrating from 0 to Infinity, so let's check a valid x value of 1000:

int2(1000, sqrt(0.245), 530, 3)
# [1] NaN

Your objective multiplies four pieces:

x <- 1000
r <- sqrt(0.245)
n <- 530
p <- 3
(1+x)^((n-1-p)/2)
# [1] Inf
(1+(1-r^2)*x)^(-(n-1)/2)
# [1] 0
x^(-3/2)
# [1] 3.162278e-05
exp(-n/(2*x))
# [1] 0.7672059

We can now see that the issue is that you're multiplying infinity by 0 (or rather something numerically equal to infinity times something numerically equal to 0), which is causing the numerical issues. Instead of calculating a*b*c*d, it will be more stable to calculate exp(log(a) + log(b) + log(c) + log(d)) (using the identity that log(a*b*c*d) = log(a)+log(b)+log(c)+log(d)). One other quick note -- the value x=0 needs a special case.

int3 = function(x, r, n, p) {
  loga <- ((n-1-p)/2) * log(1+x)
  logb <- (-(n-1)/2) * log(1+(1-r^2)*x)
  logc <- -3/2 * log(x)
  logd <- -n/(2*x)
  return(ifelse(x == 0, 0, exp(loga + logb + logc + logd)))
}
integrate(f=int3,lower=0,upper=Inf,n=530,r=sqrt(.245),p=3, stop.on.error=FALSE)
# 1.553185e+27 with absolute error < 2.6e+18

"non-finite function value" when using integrate() in R

1 Answers