R встроить возвращает неправильные результаты; настройка подразделений, rel.tol или abs.tol ничего не делает

Question

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

Я проверяю этот ответ , но моя функция во втором примере уже сосредоточена вокруг нуля.

Справочная информация. Указанные числовые интегралы соответствуют выражению для плотности суммы случайных величин и их совместной плотности is known, so we apply the formula: if Z = X + Y, then

In my model, X and Y are copula-connected, so . All distribution functions used are characterised in terms of shape and scale, and all of them have either positive or negative support, except for centred gamma (but then, the boundary is known, see below).

1. Consider the simplest example with an FGM copula (so convolution does not apply: it is for independent RV only, whereas in this case, most models do not have independent RV’s). The integration does not always work as planned. The numerical integral is often ill-behaved because the scale parameter is too small (since its estimates come from daily return data). However, when I apply a stabilising correction that seems to work for gamma distribution, centred gamma distribution begins to return nonsensical results. Consider the following example:

shapeX 

Теперь возьмем некоторые реальные значения из предыдущих оценок:

scaleX <- 0.001
scaleY <- 0.001
curve(f_Z, 0, 0.03, xlab="z", ylab="f_Z(z)", lwd=3)
integrate(f_Z, -Inf, Inf) # 0.8307993 with absolute error < 1.1e-05 --- problem!

Теперь, так как поддержка настолько мала, что квадратура не получает домен правильно, мы численно стабилизируем нашу функцию, растягивая ее, используя стабилизирующий множитель:

stab <- (scaleX+scaleY)/2 # In case one of them is close to zero, which can happen
f_Zstab <- Vectorize(function(z) integrate(function(v) f_XY(z-v*stab, v*stab), 0, Inf)$value*stab)
curve(f_Zstab, 0, 0.03, col="red", add = TRUE)
integrate(f_Zstab, -Inf, Inf) # 1 with absolute error < 3.4e-08 --- good!

Вуаля, похоже, что проблема была решена, но в других спецификациях оптимальные формы и масштабы отличаются, поэтому время от времени происходит что-то подобное.

2. Рассмотрим две переменные с центрированным гамма-распределением (то есть сглаженные; среднее значение гамма-распределенного RV - шкала Скапеса).

# Density of centred gamma with domain [shape*scale, +Inf)
dcgamma <- function(x, shape, scale) return(dgamma(x+shape*scale, shape=shape, scale=scale))
shapeX <- 25.6
shapeY <- 25.8
scaleX <- 0.007
scaleY <- 0.028
f_X <- function(x) dcgamma(x,  shape=shapeX, scale = scaleX) # Support [-shapeX*scaleX, Inf)
f_Y <- function(x) dcgamma(-x, shape=shapeY, scale = scaleY) # Support (-Inf, shapeY*scaleY]
copula <- function(x, y) 1
f_XY <- function(x, y) f_X(x)*f_Y(y)*copula(x, y)
stab <- (scaleX+scaleY)/2 # The same stabilisation technique
ulim <- shapeY*scaleY

f_Zstab <- Vectorize(function(z) integrate(function(v) f_XY(z-v*stab, v*stab), -Inf, ulim/stab)$value*stab)
curve(f_Zstab, -0.5, 0.5, xlab="z", ylab="f_Z(z)", lwd=3) # This does not look OK!
integrate(f_Zstab, -Inf, Inf) # 0.2881533 with absolute error < 9.1e-05 --- wrong!

# Now without stabilisation
stab <- 1
curve(f_Zstab, -0.5, 0.5, add=TRUE, col="red") # This looks correct!
integrate(f_Zstab, -Inf, Inf) # 0.9999983 with absolute error < 4.4e-07 --- correct!

Я решил исследовать и зафиксировать значение z=0, чтобы посмотреть на f_XY(-v*stab, v*stab) и почему его интеграл от -Inf до ulim/stab равен 0:

stab <- (scaleX+scaleY)/2
z <- 0
f0 <- function(v) f_XY(z-v*stab, v*stab)
curve(f0, -20, ulim/stab)
integrate(f0, -Inf, 0) # 70.5132 with absolute error < 0.00036
integrate(f0, -Inf, 20) # 155.4917 with absolute error < 0.0031
a <- integrate(f0, -Inf, ulim/stab) # 1.490699e-23 with absolute error < 0 --- WAT
str(a)
# List of 5
#  $ value       : num 1.49e-23
#  $ abs.error   : num 0
#  $ subdivisions: int 2
#  $ message     : chr "OK"
#  $ call        : language integrate(f = f0, lower = -Inf, upper = ulim/stab)
#  - attr(*, "class")= chr "integrate"
integrate(f0, -Inf, ulim/stab, subdivisions = 500) # Same 1.49e-23
integrate(f0, -Inf, ulim/stab, rel.tol=1e-40, abs.tol = 1e-40) # Same 1.49e-23
integrate(f0, -Inf, ulim/stab, abs.tol=0) # Same 1.49e-23

Посмотрите на эту функцию! Он хорошо себя ведет, сосредоточен вокруг 0, а его числовой интеграл по опоре все еще равен 1.49e-23!

Это очень странный результат, который вызывает много вопросов. Может быть, добавление произвольных констант поможет, но мне нужны надежные решения, поскольку в каждой функции вероятности я оцениваю 2000 таких интегралов, а затем использую стохастические максимизаторы, которые оценивают эту функцию в 2000 точках за 600 итераций. Число этих интегралов непомерно велико, и такие выбросы могут серьезно повлиять на максимизацию, поскольку в функции логарифмического правдоподобия берется log (f_Z), и 5 таких значений на 2000 приводят к отключению логарифмического правдоподобия на -250! Что можно сделать в этой ситуации? Я чувствую, что увеличение количества точек в dquagie маршрутизации Фортрана было бы решением ... но такое низкоуровневое, что его трудно реализовать, а высокоуровневые могут замедлить оценку (сотни миллионов таких интегралов).

Какое лучшее решение для стабилизации этих причудливых числовых интегралов?