Решение двойных (n-кортежных) многомерных интегралов в R

Я хотел бы написать код для решения такого рода уравнений:

уравнение для решения

Для этого я написал код ниже, однако он не решает проблему. Есть ли у вас какие-либо идеи о возможности решения такого рода интегралов в R?

t_0 = 15
mu = 0.1
lambda = 0.8
f = function(x1,x2) exp(mu*(x1+x2))*dexp(log(lambda)*(x1+x2))
f_comp = function(x2) f(x1,x2)
f_1 = function(x1) {integrate(f_comp,upper = t_0, lower = x1)}
result = integrate(f = f_1, lower = 0, upper = t_0)$value

--------- редактировать:

Учитывая ответ ниже, я адаптирую код к своему примеру, но я все же думаю, что это не правильный, по крайней мере, значение 0 для интеграла не имеет смысла.

integrate(function(x1) {
  sapply(x1, function(x1){
    integrate(function(x2)  exp(mu*(x1+x2))*dexp(log(lambda)*(x1+x2)), lower = x1, upper = t_0)$value
  })
}, 0, t_0)

кстати, я хотел бы построить общую процедуру для этого (поэтому я просто не вычисляю интеграл вручную). Это не только двойные интегралы, но и интегралы из n-кортежей, поэтому мне нужна общая процедура для такого рода вычислений.


r integral
person Francisco    schedule 26.10.2016    source источник
comment
Вы уверены, что определение функции f в R верно? У меня есть сомнения. Вы можете найти ответ в этой теме.   -  person Bhas    schedule 26.10.2016
comment
@Bhas это полезно, спасибо. Я все еще не получаю правильных значений, но, по крайней мере, в коде больше нет ошибок....   -  person Francisco    schedule 26.10.2016
comment
Я имел в виду, что, учитывая приведенную вами математическую формулу, я не думаю, что то, что вы закодировали в R как функцию f, правильно. Почему log? Где 1-....?   -  person Bhas    schedule 26.10.2016
comment
dexp определяется как [1-exp()] и exp(log(lambda)) = lambda. Однако формула не является важной частью этого вопроса.   -  person Francisco    schedule 27.10.2016
comment
Я не думаю, что dexp определяется так, как вы представляете. В справке R сказано, что dexp определяется как lambda * exp(-lambda*x)   -  person Bhas    schedule 27.10.2016
comment
ах, да, вы правы, я имею в виду pexp(). как я уже сказал, формула не важна, вычисление кратных интегралов - моя проблема. Спасибо.   -  person Francisco    schedule 27.10.2016

Ответы (1)


arrow_upward
0
arrow_downward

Нарисуйте область интеграции. Это симплекс (треугольник) с вершинами (0,0), (0,t0), (t0,t0). Для вычисления интеграла на симплексе лучше всего подойдет пакет SimplicialCubature.

t0 = 15
mu = 0.1
lambda = 0.8

library(SimplicialCubature)
f <- function(xy){
  x <- xy[1]; y <- xy[2]
  exp(-mu*(x+y)) * (1-exp(-lambda*(x+y))) 
}
S <- cbind(c(0,0), c(0,t0), c(t0,t0))
adaptIntegrateSimplex(f, S)$integral
# 29.55906

integrate(function(x1) {
  sapply(x1, function(x1){
    integrate(function(x2) exp(-mu*(x1+x2))*(1-exp(-lambda*(x1+x2))), lower = x1, 
              upper = t0)$value
  })
}, 0, t0)$value
# 29.55906
person Stéphane Laurent    schedule 03.02.2020