Вложенное интегрирование для неполной свертки плотностей Гаусса - PullRequest
Новая
       23

Вложенное интегрирование для неполной свертки плотностей Гаусса

2 голосов
/ 28 мая 2019

Пусть g(x) = 1/(2*pi) exp ( - x^2 / 2) - плотность нормального распределения со средним 0 и стандартным отклонением 1. В некоторых расчетах на бумаге появились интегралы вида

enter image description here

, где 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)

Что я получаю:

Ошибка: оценка вложена слишком глубоко: бесконечная рекурсия / Опции (выражения =)

Я понимаю, что это сложное вычисление, но для "маленького" n нечто подобное должно быть возможным.

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

1 Ответ

2 голосов
/ 28 мая 2019

Проблема заключается в том, как 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, а не ручные вычисления.

...