Triad sou.

Rのadapt package

2〜20次元関数のAdaptive Quadrature(適応的求積法)ができるパッケージ。

library(adapt)

# Example of  p - dimensional normal distribution

mu <- rbind(c(0, 0, 0))
sigma <- rbind(c(1, 0, 0), c(0, 1, 0), c(0, 0, 1))

pdfnm <- function(z) {
  1 / ( (2 * pi)^(length(z) / 2) * det(sigma)^(1/2) ) *
  exp(-0.5 * (t(z) - mu) %*% solve(sigma) %*% (z - t(mu)))
}

zlo = c(-10, -10, -10)
zup = c(0, 0, 0)
adapt(3, lo = zlo, up = zup, functn = pdfnm, eps = 1e-8)

    value       relerr     minpts     lenwrk      ifail
0.1250000 9.160987e-09      17427       3204          0

多変量標準正規分布を数値積分
時間がかかるのでlowerを-10にしてみた。