R多重積分

Question

私はこのような4つの積分持っている：この関数を計算するには

を次のように、私のRコードは次のとおりです。

tol <- 1e-05 
t.cutoff = 1e-15 

integrand <- function(kappa, sigma, eta, tau) sqrt((tau-eta)/(sigma-kappa))*exp(-(sigma-kappa)) 

integrate1 <- function(sigma, eta, tau) integrate(integrand, tol * sigma, (1 - tol) * sigma, sigma = sigma, eta = eta, tau = tau)$value 
integrate1.vec <- function(sigma, eta, tau) mapply(integrate1, sigma, eta, tau) 

integrate2 <- function(eta, tau) integrate(integrate1.vec, tol * eta, (1 - tol) * eta, eta = eta, tau = tau)$value 
integrate2.vec <- function(eta, tau) mapply(integrate2, eta, tau) 

integrate3 <- function(tau) integrate(integrate2.vec, tol * tau, (1 - tol) * tau, tau = tau)$value 
integrate3.vec <- function(tau) mapply(integrate3, tau) 


c_t <- function(t) ifelse(t > t.cutoff, integrate(integrate3.vec, tol * t, (1 - tol) * t)$value, 0) 

c_t(1) 

私は上記のコードは、効率的な低だと思います;この種のインテグラルには改善や利用可能なパッケージがあるのだろうかと思います。前もって感謝します！

Answer 1

この統合領域は、cumsumをシンプレックスに適用することで取得できます。したがって、SimplicialCubatureパッケージを使用することができます。

library(SimplicialCubature) 
integrand <- function(y){ 
    x <- cumsum(y) 
    sqrt((x[4]-x[3])/(x[2]-x[1]))*exp(-(x[2]-x[1])) 
} 
t <- 1 
S <- t*CanonicalSimplex(4) 
adaptIntegrateSimplex(integrand, S, maxEvals = 100000L, absError = 1e-4) 

結果：

> adaptIntegrateSimplex(integrand, S, maxEvals = 100000L, absError = 1e-4) 
$integral 
[1] 0.05961484 

$estAbsError 
[1] 9.991182e-05 

私はMathematicaのでチェックしたのだが、0.0596353を与えます。

我々はabsErrorを減らすことによって、この結果に近づくが、これはmaxEvalsを増やす必要があります。

> adaptIntegrateSimplex(integrand, S, maxEvals = 1000000L, absError = 1e-5) 
$integral 
[1] 0.05963446 

$estAbsError 
[1] 9.993167e-06