4D Монте-Карло Интеграция R - PullRequest
1 голос
/ 02 февраля 2020

Как интегрировать функцию, используя метод Монте-Карло

F(X[1],X[2],X[3],X[4])

Что зависит от 4 переменных по 4 измерениям?

Я имею в виду int_{0}^{1} int_{0}^{1} int_{0}^{1} int_{0}^{1} X[1]X[2]X[3]X[4] dX[1] dX[2] dX[3] dX[4]

Функция UPD

data1 =  rnorm(100, 0, 1)
data2 =  rnorm(100, 0.1, 0.5)
data3 =  rnorm(100, 0.2, 0.8)
data4 =  rnorm(100, 0.3, 0.9)

kernel1 = kdensity(data1,kernel = 'gaussian')
kernel2 = kdensity(data2,kernel = 'gaussian')
kernel3 = kdensity(data3,kernel = 'gaussian')
kernel4 = kdensity(data4,kernel = 'gaussian')


f <- function(X) {
  return(X[1]*kernel1(X[1])*kernel2(X[2])*kernel3(X[3])*kernel4(X[4]))
}

и я хочу интегрировать его

int_{0}^{1} int_{0}^{1} int_{0}^{1} int_{0}^{1} f dX[1] dX[2] dX[3] dX[4]

1 Ответ

1 голос
/ 03 февраля 2020

Благодаря вашей особой структуре интегрального выражения вы можете переписать вложенный интеграл в произведение интегралов.

Таким образом, приведенное ниже решение может быть примером (со случайным начальным числом set.seed(1)) для вас:

g <- Vectorize(function(ker) {integrate(ker,0,1)}$value)
gE <- function(ker) {integrate(function(x) x*ker(x),0,1)}$value
res1 <- prod(c(gE(kernel1),g(c(kernel2,kernel3,kernel4))))

таким, что

> res1
[1] 0.01559343

Вложенный Integral : вам нужно сначала переписать вашу функцию f, то есть

f <- function(x1,x2,x3,x4) {
  return(x1*kernel1(x1)*kernel2(x2)*kernel3(x3)*kernel4(x4))
}

res2 <- integrate(Vectorize(function(x4) 
  integrate(Vectorize(function(x3,x4) 
    integrate(Vectorize(function(x2,x3,x4) 
      integrate(f,0,1,x2,x3,x4)$value),
      0,
      1,
      x3,
      x4)$value),
    0,
    1,
    x4)$value),
  0,
  1)$value

так, чтобы

> res2
[1] 0.01559343
...