"Non-Finite Function Value" When Using Integrate() in R

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

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 error in 'integrate' function in R

The problem is that the function

exp(y)/(1+exp(y)) 

is rounded to NaN when y is too big. You can avoid this replacing it with 1 when y is too big. This function will play the trick:

fun5<-function(y,mu=mu0,lsig=-lsig0) {
res = ifelse(y<100, exp(y)/(1+exp(y)) * 1/sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2),
             1/sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2))
return(res)}

and now this will work

integrate(fun5,-Inf,Inf,mu=2.198216,lsig=-3)$value
[1] 0.9

Non-finite function value with integrate() R although solution exists

In line with @Bhas's answer, I would go for the following solution:

> f <- function(x){ifelse(x!=1,((1.-x)^2)/(abs(1.-x))^(1/3),0)} # Set f(1)=0 since it is the limit of 'f' at 1.
> integrate(f,lower=0,upper=1.6,abs.tol=1E-7)
0.4710361 with absolute error < 2.2e-08

'ifelse' avoids problems related to a vectorized 'x'

R integrate() function returns 'non-finite function value'

You should be checking interpolated values of your function, not the integration domain. One of the values must be NaN, and I can't provide more information because there is no data.

I can only guess what is going on. Since you are using approxfux, the interpolated function (latDensFunc) returns NaN for one of the points. This is what leads to NaNs in lonSq.

Numerical Solution for Integral between -Inf and Inf: Error non-finite function value

There are a couple of things to point out here. Firstly, your function needs to have as its first argument the variable over which you want to integrate, so you need to rewrite your function as:

random_walk_func<-function(x, t, A, sigma, y)
{
  a1 <- (2*A/(sigma))*exp((4*A*(y-x+(4*A*t)))/(sigma))
  b1 <- erfc((y-x+(8*A*t))/(2*sqrt(sigma*t)))
  a1 * b1
}

Secondly, remember that this is numeric rather than symbolic integration, so you need to have values for all the other parameters you are passing to your function. I have no idea what you want these to be, so let's set them all to 1:

t <- A <- sigma <- y <- 1

Thirdly, it's a good idea to look at what you're integrating if your are getting infinity errors. If there are infinite values among the evaluated points, then you will get an error rather than a numeric result:

x <- seq(-10, 10, 0.01)
 
plot(x, random_walk_func(x, t, A, sigma, y), type = "l")

We can see that we will get an excellent approximation of the integral if we choose limits of -10 and 10:

integrate(random_walk_func, lower = -10 , upper = 10, 
          t = t, A = A, sigma = sigma, y = y)$value
#> [1] 1

However, ultimately the reason why you are getting the error is that a1 gets monstrously large very quickly the further from the central peak that we go, and b1 becomes infintesimal. Even though their product is nearly zero, the intermediate calculations are beyond R's numerical tolerance, which is what breaks the calculation. Once a1 exceed about 10^308, R will call it Inf and a1 * b1 is therefore also Inf.

The way round this is to calculate a1 and b1 as logs, then return their exponentiated sum. So if you do:

random_walk_func <- function(x, t, A, sigma, y)
{
  a1 = log(2 * A / sigma) + 4 * A * (y - x + (4 * A * t)) / sigma
  b1 = log(erfc((y - x + 8 * A * t) / (2 * sqrt(sigma * t))))
  exp(a1 + b1)
}

Then you get:

integrate(random_walk_func, lower = -Inf, upper = Inf, 
          t = t, A = A, sigma = sigma, y = y)$value
#> [1] 1

How to solve this non-finite function value in R?

I gave it a try - the problem is that integrate expects the handed over function to be able to deal with input vectors and output a vector of same size.

Luckily the solution is easy - just wrap your function in sapply.

The following code works:

f <- function(tau) {integrate(function(tau1)sqrt(1/(tau-tau1)),lower=0.01*tau,upper=0.99*tau)$value}

intfun <- function(x) sapply(x,f) 

C1_B <- function(T){integrate(function(tau) intfun(tau),lower=0.01*T,upper=0.99*T)$value}
C1_B(0.5)

Integration in R giving non-finite function value

I would approach it using the polynom package, which can calculate the derivative for us, and then apply the formula from Wolfram.

revol_coef <- 2.393794315 * c(-2953570085669 / 16500,
                              835541173981 / (8250 * 2.24581),
                              -7563546913/ (330 * 2.24581^2),
                              171163181 / (66 * 2.24581^3),
                              -1613880 / (11 * 2.24581^4),
                              3320 / 2.24581^5)
y <- polynomial(revol_coef)
y_d <- deriv(y)

f <- function(x) {
    2 * pi * predict(y, x) * sqrt(1 + predict(y_d, x)^2)
}

integrate(f, lower = 19.538547, upper = 20.3245805)

The output I get is

-45.71118 with absolute error < 0.0017

The answer is negative, due to the surface being specified in the lower-right quadrant of the x-y plane (as per your plot). To fix this, take the negative of the number supplied to get a surface area.

Translating formulas from Wolfram Alpha, manipulating them, and then putting them into R can easily involve making a mistake.

In this case, I would recommend doing as much of the work within R, as the tools/functions to do the necessary calculations and simplifications are readily available.

---

Related Topics