Интегрируйте гладкую функцию, когда интегрировать () не удается
При попытке интегрировать следующую гладкую функцию с integrate() в R я всегда получаю сообщение об ошибке: либо «интеграл, вероятно, расходится», либо «достигнуто максимальное количество подразделений»:
Проблема, по-видимому, в том, что функция приближается к нулю примерно на +/-1, но трудно правильно установить пределы, и увеличение subdivisions или уменьшение rel.tol не помогает.
Существуют ли другие процедуры интеграции в R (или в пакетах R), которые я мог бы попробовать? Есть ли, например. подпрограмма для простой интеграции Симпсона, доступная для функций с равноудаленной выборкой?
Обновлено: как кто-то спросил, вот реализация функции:
a <- -1
b <- 1
sigma <- 0.1
# P(Y >= t)
PY <- function(t) {
za <- (t-a)/sigma
zb <- (t-b)/sigma
1 + sigma/(b-a) * (zb*pnorm(zb) + dnorm(zb) - za*pnorm(za) - dnorm(za))
}
# P(Y >= t | X=x)
PY.x <- function(t, x) {
1 - pnorm((t-x)/sigma)
}
# density of Y: p(y)
p.y <- function(y) {
zb <- (b-y)/sigma
za <- (a-y)/sigma
1/(b-a) * (pnorm(zb) - pnorm(za))
}
# Var in numerator
Var.Yt <- function(t) {
integrand <- function(x) {
(PY.x(t,x) - PY(t))^2
}
integrate(integrand, a, b)$value
}
# THIS IS THE FUNCTION TO BE INTEGRATED
f <- function(t) {
Var.Yt(t)*p.y(t)
}
Вычисление Var.Yt(t) выполняется успешно для любого значения t, но интегрирование f завершается с ошибкой, упомянутой выше.
Это немного сложно, потому что сама функция определяется вызовом integrate, однако я проверил, что эта внутренняя интеграция всегда работает для произвольных параметров. Я добавлю это к моему вопросу.
Вы можете попробовать pracma::integral()
(1) Зависит от того, над какой целью вы работаете, но в любом случае, учитывая определение в терминах гауссовой плотности и ее интеграла (который пропорционален erf), может быть возможно разработать символическое выражение для рассматриваемого интеграла. по сигме. (2) Похоже, вы вычисляете условную дисперсию для распределения Гаусса. Если это так, я думаю, что это имеет известный результат. См., например. Ричард фон Мизес, «Математическая теория вероятностей и статистика», глава VIII, если я правильно помню. Несомненно, такие результаты представлены и в других местах.
@robert-dodier Да, я думаю, что интеграл можно вычислить в символической форме, потому что подынтегральное выражение представляет собой многочлен от x, phi(x) и Phi(x) с членом высшего порядка x^2 phi(x)^2 Phi(x)^2, чтобы его можно было вычислить путем частичного интегрирования и формул из en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions. В настоящее время у меня не было времени на это, поэтому я прибегнул к «быстрому пути» численного интегрирования ценой некоторой потери точности. Преимущество численного интегрирования в том, что я могу адаптировать код для более сложных вариантов использования.
Ответы 1
Я думаю, вы забыли векторизовать свою функцию f, прежде чем отправить ее в integrate или pramca::integral.
Здесь вы можете увидеть, что они будут работать, если мы подадим заявку Vecotorize()
> pracma::integral(Vectorize(f), -2, 2)
[1] 0.2799594
> integrate(Vectorize(f), -2, 2)
0.2799594 with absolute error < 5.4e-06
Ах, да, вот оно! Теперь я даже могу использовать лимиты -Inf и Inf. Спасибо!
показать свою гладкую функцию