Coding a multiple integral function in R

问题:

With the goal of turning the following into a function, I was wondering how I can write the following double integral in terms of R codes?: ($\bar{x} = \mu$):

回答1:

Assuming pi0 and pi1 implement your functions $\pi_0$ and $\pi_1$ in a vectorized way, a possible solution is:

integral <- function(n, mu, s, pi0, pi1) {

  C <- (2 * pi)^(-n/2) 
  C * integrate(f = function(sigmavec) sapply(sigmavec, function(sigma) {
    integrate(f = function(delta) {

      exp(-n/2 * ((mu / sigma - delta)^2 + (s / sigma)^2)) * pi1(delta)

    }, lower = -Inf, upper = Inf)$value

  }) * pi0(sigmavec) / (sigmavec^n), lower = 0, upper = Inf)$value

}

# Tests
integral(n = 1, mu = 0, s = 1, pi0 = dnorm, pi1 = dnorm)
# [1] 0.0473819
integral(n = 1, mu = 0, s = 1, pi0 = function(sigma) 1/sigma, pi1 = dcauchy)
# [1] 0.2615783

回答2:

Note sure if this question is on topic, but I am open to answer.

May be you should ask a more general question, how to write/computing integral using computer program (code)? There at least are two ways

1. Using numerical integration, such as Monte Carlo method
2. Using symbolic toolbox to solve the problem analytically and plugin the numerical value.

Examples on $\int_0^1 x^2$

f<-function(x){
  x^2
}
curve(f,0,1)

# method 1
integrate(f,lower=0,upper = 1)

# method 2
library(Ryacas)
x <- Sym("x")
f <- function(x) {
  x^2
}
f2=yacas(yacas(Integrate(f(x), x)))
f2

x <- 1
Eval(f2)