目标是通过这些随机变量之一的密度函数来计算 n 个 IID 随机变量之和的密度函数:
以下是我对此的天真的尝试:
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 和产生适当缩放和适当移位的最终密度函数?
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)
)