Вложенная интеграция для неполной свертки гауссовых плотностей
Пусть g(x) = 1/(2*pi) exp ( - x^2 / 2) будет плотностью нормального распределения со средним значением 0 и стандартным отклонением 1. В некоторых расчетах на бумаге появились интегралы вида
где c>0 — положительное число.
Поскольку я не мог оценить это вручную, у меня возникла идея приблизить и построить его. Я попробовал это в R, потому что R предоставляет функцию dnorm и функцию для вычисления интегралов.
Вы видите, что мне нужно численно интегрировать n раз, где n должно быть выбрано вызовом функции построения графика. В моем коде есть цикл for для итеративного создания этих «неполных» сверток.
Например, даже при n=3 и c=1 это дает мне ошибку. n=2 (таким образом, это одна интеграция) работает.
N = 3
ngauss <- function(x) dnorm(x , mean = 0, sd = 1)
convoluts <- list()
convoluts[[1]] <- ngauss
for (i in 2:N) {
h <- function(y) {
g <- function(z) {ngauss(y-z)*convoluts[[i-1]](z)}
return(integrate(g, lower = -1, upper = 1)$value)
}
h <- Vectorize(h)
convoluts[[i]] <- h
}
convoluts[[3]](0)
Что я получаю:
Error: evaluation nested too deeply: infinite recursion / options(expressions=)?
Я понимаю, что это сложное вычисление, но для «маленьких» n нечто подобное должно быть возможно.
Может быть, кто-то может помочь мне исправить мой код или дать рекомендацию, как я могу реализовать это лучше. Другой язык, который более подходит для этого, также будет в порядке.
Ответы 1
Проблема, по-видимому, заключается в том, как integrate работает с переменными в разных средах. В частности, он не совсем правильно обрабатывает i на каждой итерации. Вместо этого используя
h <- evalq(function(y) {
g <- function(z) {ngauss(y - z) * convoluts[[i - 1]](z)}
integrate(g, lower = -1, upper = 1)$value
}, list(i = i))
делает свое дело и, скажем, настройка N <- 6 быстро дает
convoluts[[N]](0)
# [1] 0.03423872
Поскольку ваше интегрирование — это просто PDF суммы N независимых стандартных нормалей (за которой следует N(0, N)), мы также можем проверить этот подход, установив lower = -Inf и upper = Inf. Тогда с N <- 4 мы имеем
dnorm(0, sd = sqrt(N))
# [1] 0.1994711
convoluts[[N]](0)
# [1] 0.1994711
Итак, для практических целей, когда c = Inf, вам лучше использовать dnorm, чем ручные вычисления.
Честно говоря, я думаю, что теперь это вычисление мне не поможет. Даже для «N<-6» это занимает слишком много времени
@Falrach, возможно, вы уже изучали это, но для общего c у нас есть сумма симметрично усеченных нормальных случайных величин, поэтому, возможно, jstor.org/стабильный/2236545 или какая-либо другая статья могут дать хорошее приближение.
на самом деле я не смотрел на это. Я думаю, что эти понятия связаны с моей проблемой, но мой интеграл больше соответствует случаю усечения после каждой свертки, чем усечению до свертки, как в этой статье. Но я думаю, что понятие усеченных случайных величин поможет мне выяснить, есть ли уже что-то, связанное с моим вопросом.
Прежде всего, спасибо за ответ. Я реализовал ваш код, и, по крайней мере, больше не получаю ошибок. Кажется, это работает. Конечно, меня интересует не только оценка свертки[[i]] на нуле, но и построение ее по всей реальной линии. Если c = Inf, конечно, я мог бы лучше использовать полугрупповое свойство нормального распределения в отношении свертки. Путем повторной свертки масса вытекает из каждого компакта. Но решающим моментом здесь является то, что мы не интегрируем свертку по всей реальной прямой. Может произойти изменение поведения в связи с этой потерей концентрации.