Я получил очень маловероятную, но очень опасную числовую ошибку при объединении тысяч функций с хорошим поведением в R с помощью встроенной функции integrate.
История (можно пропустить). Моя проблема связана с максимальной вероятностью и основана на сильно нелинейной функции (из 10–20 параметров), для которой аналитическое выражение не существует, поэтому это требует вычисления тысяч интегралов для одной оценки. Я создал MWE, который содержал эту ошибку. Для оптимизации этой функции из-за нескольких локальных оптимумов я пробую 1000 точек на 1000 итераций (с помощью методов без производных, таких как рой частиц от hydroPSO и дифференциальное развитие от DEoptim), поэтому для одной модели я приходится вычислять более миллиарда интегралов (!), и существует 200 моделей-кандидатов, каждая из которых требует более поздней переоценки, поэтому общее число интегралов превышает триллион. Я хотел бы найти самое быстрое решение, которое дает достаточную точность.
Функция является произведением двух функций плотности (гамма или подобных), умноженных на некоторое положительное выражение, и плотность суставов вычисляется по формуле f_{X+Y}(z) = int_{supp Y} f_{X+Y}(z-y, y) dy. Я не могу использовать свертку, потому что X и Y не являются независимыми в общем случае. Поддержка Y в моем случае - (-Inf, 0]. Параметр масштаба распределения очень мал (модель похожа на GARCH), поэтому очень часто стандартная процедура интеграции не может интегрировать ненулевую функцию в очень маленький участок отрицательной линии (например, * 1017). * где он принимает огромные значения и 0 везде, где он очень старается вычислить квадратуру), а R integrate часто возвращает машинный эпсилон, потому что он не может найти точки в этом диапазоне, где функция принимает значения намного большие чем ноль. Чтобы бороться с этой проблемой, я растягиваю функцию вокруг нуля с помощью обратного параметра масштаба, вычисляю интеграл и затем делю его на масштаб, т.е. е. integrate(f(x/scale)/scale$value). Однако иногда такое повторное масштабирование также не удавалось, поэтому я реализовал проверку безопасности, чтобы увидеть, является ли значение масштабированной функции подозрительно низким (т. Е. <1e-8), и затем заново вычислил интеграл. Масштабируемая функция работала как талисман, возвращая хорошие значения, когда немасштабируемая функция не работала, а в редких случаях функция с измененным масштабом возвращала машинный эпсилон, а не масштабируемая работала.
До сегодняшнего дня, когда интеграция масштабированной функции неожиданно дала значение 1,5 вместо 3,5. Конечно, функция прошла проверку безопасности (потому что это правдоподобное значение, а не машинный эпсилон и некоторые другие значения были меньше этого, поэтому было в общем диапазоне). Оказалось, примерно в 0,1% всех случаев integrate недооценили функцию. MWE ниже.
Сначала мы определим функцию x и необязательный параметр numstab, который определяет масштабирование.
cons <- -0.020374721416129591
sc <- 0.00271245601724757383
sh <- 5.704
f <- function(x, numstab = 1) dgamma(cons - x * numstab, shape = sh, scale = sc) * dgamma(-x * numstab, shape = sh, scale = sc) * numstab
Далее мы строим график, чтобы убедиться, что масштабирование работает правильно.
curve(f, -0.06, 0, n = 501, main = "Unscaled f", bty = "n")
curve(f(x, sc), -0.06 / sc, 0, n = 501, main = "Scaled f", bty = "n")
И затем мы проверяем этот интеграл суммированием:
sum(f(seq(-0.08, 0, 1e-6))) * 1e-6 # True value, 3.575294
sum(f(seq(-30, 0, 1e-4), numstab = sc)) * 1e-4 # True value, 3.575294
str(integrate(f, -Inf, 0)) # Gives 3.575294
# $ value : num 3.58
# $ abs.error : num 1.71e-06
# $ subdivisions: int 10
str(integrate(f, -Inf, 0, numstab = sc))
# $ value : num 1.5 # WTF?!
# $ abs.error : num 0.000145 # WTF?!
# $ subdivisions: int 2
Это остановится только в двух подразделениях! Теперь, чтобы увидеть, что происходит во время интеграции, мы создаем глобальный объект и обновляем его каждый раз, когда процедура интеграции что-то делает.
global.eval.f <- list()
f.trace <- function(x, numstab = 1) {
this.f <- f(x, numstab)
global.eval.f[[length(global.eval.f) + 1]] <<- list(x = x, f = this.f)
return(this.f)
}
integrate(f.trace, -Inf, 0)
Теперь мы визуализируем этот процесс интеграции.
library(animation)
l <- length(global.eval.f)
mycols <- rainbow(l, end = 0.72, v = 0.8)
saveGIF({
for (i in 1:l) {
par(mar = c(4, 4, 2, 0.3))
plot(xgrid <- seq(-0.1, -0.01, length.out = 301), f(xgrid), type = "l", bty = "n", xlab = "x", ylab = "f(x)", main = "Function without stabilisation")
for (j in 1:(l2 <- length(this.x <- global.eval.f[[i]]$x))) lines(rep(this.x[j], 2), c(0, global.eval.f[[i]]$f[j]), col = mycols[i], type = "b", pch = 16, cex = 0.6)
legend("topleft", paste0("Quadrature: ", i), bty = "n")
text(rep(-0.1, l2), seq(325, 25, length.out = l2), labels = formatC(sort(this.x), format = "e", digits = 2), adj = 0, col = ifelse(sort(this.x) > -0.1 & sort(this.x) < -0.01, mycols[i], "black"), cex = 0.9)
}
}, movie.name = "stab-off-quad.gif", interval = 1 / 3, ani.width = 400, ani.height = 300)
И то же самое, но в другом масштабе.
global.eval.f <- list()
integrate(f.trace, -Inf, 0, numstab = sc)
l <- length(global.eval.f)
mycols <- rainbow(l, end = 0.7, v = 0.8)
saveGIF({
for (i in 1:l) {
par(mar = c(4, 4, 2, 0.3))
plot(xgrid <- seq(-0.1 / sc, -0.01 / sc, length.out = 301), f(xgrid, sc), type = "l", bty = "n", xlab = "x", ylab = "f(x)", main = "Function with stabilisation")
for (j in 1:(l2 <- length(this.x <- global.eval.f[[i]]$x))) lines(rep(this.x[j], 2), c(0, global.eval.f[[i]]$f[j]), col = mycols[i], type = "b", pch = 16, cex = 0.6)
legend("topleft", paste0("Quadrature: ", i), bty = "n")
text(rep(-0.1 / sc, l2), seq(325 * sc, 25 * sc, length.out = l2), labels = formatC(sort(this.x), format = "e", digits = 2), adj = 0, col = ifelse(sort(this.x) > -0.1 / sc & sort(this.x) < -0.01 / sc, mycols[i], "black"), cex = 0.9)
}
}, movie.name = "stab-on-quad.gif", interval = 1 / 3, ani.width = 400, ani.height = 300)
Проблема в том, что я не могу попробовать различные стабилизирующие множители для функции, потому что мне приходится вычислять этот интеграл триллион раз, поэтому даже в кластере суперкомпьютеров это занимает недели. Кроме того, уменьшение rel.tol до 1e-5 немного помогло, но я не уверен, гарантирует ли это успех (а уменьшение до 1e-7 в некоторых случаях замедляло вычисления). И я посмотрел на квадратурный код Фортрана, чтобы увидеть правило интеграции.
Время можно увидеть ниже (я добавил дополнительную попытку с меньшим допуском).
Как я могу убедиться, что процедура интеграции не даст таких неправильных результатов для такой функции, и интеграция все еще будет быстрой?