Search code examples
rnumerical-integrationintegral

How to calculate integral inside an integral in R?

I need to evaluate an integral in the following form:

\int_a^b f(x) \int_0^x g(t)(x-t)dtdx

Can you please suggest a way? I assume that this integral can't be done in the standard approach suggested in the following answer:

Standard approach

Update: Functions are added in the following image. f(x) basically represents a pdf of a uniform distribution but the g(t) is a bit more complicated. a and b can be any positive real numbers.

enter image description here

Solution

  • The domain of integration is a simplex (triangle) with vertices (a,a), (a,b) and (b,b). Use the SimplicialCubature package:

    library(SimplicialCubature)

alpha <- 3
beta <- 4
g <- function(t){
  ((beta/t)^(1/2) + (beta/t)^(3/2)) * exp(-(t/beta + beta/t - 2)/(2*alpha^2)) / 
    (2*alpha*beta*sqrt(2*pi))
}
a <- 1
b <- 2
h <- function(tx){
  t <- tx[1]
  x <- tx[2]
  g(t) * (x-t)
}

S <- cbind(c(a, a), c(a ,b), c(b, b))
adaptIntegrateSimplex(h, S)
# $integral
# [1] 0.01962547
# 
# $estAbsError
# [1] 3.523222e-08

    Another way, less efficient and less reliable, is:

    InnerFunc <- function(t, x) { g(t) * (x - t) }
InnerIntegral <- Vectorize(function(x) { integrate(InnerFunc, a, x, x = x)$value})
integrate(InnerIntegral, a, b)
# 0.01962547 with absolute error < 2.2e-16