Search code examples
rmathnumerical-methodsintegral

Integrate a smooth function when integrate() fails

When trying to integrate the following smooth function with integrate() in R, I always obtain an error, either the "integral is probably divergent" or "maximum number of subdivisions reached":

enter image description here

The problem seems to be that the function approaches zero at around +/-1, but it is difficult to set the limits appropriately, and increasing subdivisions or reducing rel.tol does not help.

Are there other integration routines in R (or in R packages) that I could try? Is there e.g. a routine for simple Simpson integration available for equidistantly sampled functions?

Edit: As someone asked, here is the implementation of the function:

a <- -1
b <- 1
sigma <- 0.1

# P(Y >= t)
PY <- function(t) {
  za <- (t-a)/sigma
  zb <- (t-b)/sigma
  1 + sigma/(b-a) * (zb*pnorm(zb) + dnorm(zb) - za*pnorm(za) - dnorm(za))
}

# P(Y >= t | X=x)
PY.x <- function(t, x) {
  1 - pnorm((t-x)/sigma)
}

# density of Y: p(y)
p.y <- function(y) {
  zb <- (b-y)/sigma
  za <- (a-y)/sigma
  1/(b-a) * (pnorm(zb) - pnorm(za))
}

# Var in numerator
Var.Yt <- function(t) {
  integrand <- function(x) {
    (PY.x(t,x) - PY(t))^2
  }
  integrate(integrand, a, b)$value
}

# THIS IS THE FUNCTION TO BE INTEGRATED
f <- function(t) {
  Var.Yt(t)*p.y(t)
}

The computation of Var.Yt(t) succeeds for any value of t, but the integration of f fails with the aforementioned error.

Solution

  • I think you forgot to vectorize your function f before pushing it to integrate or pramca::integral.

    Here you can see they would work if we apply Vectorize()

    > pracma::integral(Vectorize(f), -2, 2)
[1] 0.2799594

> integrate(Vectorize(f), -2, 2)
0.2799594 with absolute error < 5.4e-06