高斯密度不完全卷积的嵌套积分

Question

设g(x) = 1/(2*pi) exp ( - x^2 / 2)为均值为0、标准差为1的正态分布的密度。在纸上的一些计算中出现了

形式的积分

其中 c>0 为正数。

由于我无法用手评估它，所以我有了近似并绘制它的想法。我在 R 中试过这个，因为 R 提供了 dnorm 函数和一个积分函数。

你看我需要对n次进行数值积分，其中n是通过调用plot函数来选择的。我的代码有一个 for 循环来迭代地创建那些 "incomplete" 卷积。

例如，即使 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=)?

我知道这是一项艰巨的计算，但对于 "small" n 来说，类似的事情应该是可能的。

也许有人可以帮助我修复我的代码或提供建议以更好地实现它。另一种更适合这个的语言也可以。

Answer 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

由于您的积分只是 N 独立标准法线总和的 pdf（然后遵循 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 而不是手动计算。

高斯密度不完全卷积的嵌套积分

Nested integration for incomplete convolution of gauss densities

nested

r

convolution