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