R: double integral with dependencies

Question

My goal is to compute the following double integral in R.

I looked at the previous solutions, such as double integral in R. Following the solution by G5W, I came up with the code

inner_func <- function(x) { 
  alpha=23
  beta=14
  return(x^(alpha-1)*(1-t-x)^(beta-1))
}

innerintegral <- Vectorize(
  function(t) {
    integrate(inner_func,0,1-t)$value
  }
)

integrate(innerintegral,0,1)

This does not work. I think because the inner function itself is dependent on the limit, I am not getting any solution.

I also looked into solution by MrFlick and ran the following code, which does give me an output.

fun0 <-  function(x,t){
  alpha <- 10
  beta <- 10
  return(x^(alpha-1)*(1-t-x)^(beta-1))
}

  integrate(function(t) { 
    sapply(t, function(t) {
      integrate(function(x) fun0(x,t), 0, 1-t)$value
    })
  }, 0, 1)$value

[1] 5.412544e-08

I am not sure if this is the right way to do it or even the solution is correct. Please let me know if this is the right procedure and the solution is correct or not.

Answer 1

@MrFlick's method is correct. I find the same result with the SimplicialCubature package:

library(SimplicialCubature)

f <- function(xt) {
  x <- xt[1]; t <- xt[2]
  alpha <- 10
  beta  <- 10
  x^(alpha-1) * (1-t-x)^(beta-1)
}

# simplex
S <- cbind(c(0, 0), c(1, 0), c(0, 1))

adaptIntegrateSimplex(f, S)
# $integral
# [1] 5.412543e-08
# 
# $estAbsError
# [1] 5.071444e-13

I think the simplicial method is better, because the method you use ignores the errors on the inner integral.

Bonus

Here is the value of the integral without integrating:

> beta(10,10) / 20
[1] 5.412544e-08

The general formula is beta(alpha, beta) / (alpha + beta).