在 R 中使用 FFT 确定 IID 和的密度函数

Question

目标是通过以下随机变量之一的密度函数计算 n 个 IID 随机变量之和的密度函数：

1. 通过 fft 将密度函数转换为特征函数
2. 将特征函数提升到 n
3. 变换所得特征函数通过 fft(inverse=TRUE) 进入感兴趣的密度函数

下面是我对此的天真的尝试：

sum_of_n <- function(density, n, xstart, xend, power_of_2)
{
  
  x <- seq(from=xstart, to=xend, by=(xend-xstart)/(2^power_of_2-1))
  y <- density(x)
  fft_y <- fft(y)
  fft_sum_of_y <- (fft_y ^ n) 
  sum_of_y <- Re(fft(fft_sum_of_y, inverse=TRUE))
  return(sum_of_y)
  
}

在上面，密度是任意密度函数：例如

density <- function(x){return(dgamma(x = x, shape = 2, rate = 1))} 

n 表示被求和的 IID 随机变量的数量。 xstart 和 xend 是随机变量的近似支持的开始和结束。 power_of_2 是所使用的数值向量的 2 长度的幂。据我了解，二次幂的长度提高了 fft 算法的效率。

我至少部分理解为什么上述内容通常不能按预期工作。首先，值本身将无法正确缩放，因为 fft(inverse=TRUE) 默认情况下不会标准化。然而，我发现当我除以向量的长度时，这些值仍然不正确，即

sum_of_y <- sum_of_y / length(sum_of_y)

基于我对 fft 的有限理解是归一化计算。其次，由于执行 fft 时发生的零频率偏移（如果我错了，请有人纠正我），所得矢量将异相。例如，我尝试使用 pracma 的 fftshift 和 ifftshift，但它们似乎无法正确解决此问题。对于对称分布（例如正态分布），这并不难解决，因为相移通常恰好是一半，因此类似的操作

sum_of_y <- c(sum_of_y[(length(y)/2+1):length(y)], sum_of_y[1:(length(y)/2)])

可以作为校正。然而，对于像上面的伽玛分布这样的不对称分布，这种方法会失败。

总之，是否对上述代码进行了调整，从而为 IID 和产生适当缩放和适当移动的最终密度函数？

Answer 1

dsumf <- function(d, n, a, b, k) {
  x <- seq(a, b, length.out = k)
  p <- d(x)
  s <- sum(p)
  data.frame(x = x, d = Re(fft(fft(p/s)^n, TRUE))*s/k)
}

# use fft to approximate the density of `sum(rgamma(5, 2, 1))`
dsum <- dsumf(\(x) dgamma(x, 2, 1), 5, 0, 32, 2^10)
# plot the approximation against the true density
plot(dsum$x, dgamma(dsum$x, 10, 1), type = "l", col = "blue",
     xlab = "x", ylab = "density")
lines(dsum$x, dsum$d, lty = 2, col = "orange")
legend(
  "topright",
  legend = c("True density", "FFT approximation"),
  col = c("blue", "orange"),
  lty = c(1, 2)
)