使用FFT对独立随机变量求和

Question

基于对另一个问题的非常有用的答案，我编写了一个函数来计算随机变量之和的密度函数。该函数采用两个密度函数

f1

和

f2

，然后分别生成具有这些密度的

n1

和

n2

随机变量之和的 PDF 之和（因此对

n1+n2

独立随机变量求和）：

dsumf2 <- function(f1, f2, n1=1, n2=1, a, b, k=2^14) {
  # Perform FFT
  x <- seq(a, b, length.out = k)
  p1 <- f1(x)
  p2 <- f2(x)
  s1 <- sum(p1)
  s2 <- sum(p2)
  d0 = fft(fft(p1/s1)^n1*fft(p2/s2)^n2, TRUE)
  d0 = Re(d0) * exp(n1/(n1+n2)*log(s1) + n2/(n1+n2)*log(s2)-log(k))
  data.frame(
    x = x, 
    d = d0
  )
}

当最小值为 0 时，该函数似乎运行良好：

## Works well; test with Monte Carlo
## Sum 2 gamma(1,5) and 1 gamma(2, 10)
x = rgamma(100000,shape = 1, scale = 5) + rgamma(100000,shape = 1, scale = 5) + rgamma(100000,shape = 2, scale = 10)
hist(x, freq=FALSE, breaks = 100)
dsum <- dsumf2(\(x) dgamma(x, shape = 1, scale = 5), \(x) dgamma(x, shape = 2, scale = 10), n1=2, a=0, b=150)
lines(dsum$x, dsum$d, col = 'red')

但是，当最小值不为 0 时，它会失败，如本例所示：

x = rnorm(100000) + rnorm(100000) + rnorm(100000,-5,5)
hist(x, freq=FALSE, breaks = 100)
dsum <- dsumf2(\(x) dnorm(x), \(x) dnorm(x, -5, 5), n1 = 2, a = min(x)-10, b = max(x)+1)
lines(dsum$x, dsum$d, col = 'red')

生成的 PDF 环绕，这与 this Question 中的问题类似，所以我假设这是同一个问题。

plot(dsum$x, dsum$d, col = 'red', typ='l')

我已经尝试使用

SynchWave::ifftshift

和

SynchWave::fftshift

，如该答案所示，但我无法让它正常工作。有没有办法编写函数来以正确的方式转换 PDF，而不管底层发行版及其支持如何？

Answer 1

我看到了一些非常有帮助的评论，但没有回复来知道您是否完成了这项工作。但是，为了解决函数 dsumf2 中的周期性问题，您确实可以利用 fftshift 操作来正确对齐傅立叶变换密度。

dsumf2 <- function(f1, f2, n1=1, n2=1, a, b, k=2^14) {
  # Perform FFT
  x <- seq(a, b, length.out = k)
  p1 <- f1(x)
  p2 <- f2(x)
  s1 <- sum(p1)
  s2 <- sum(p2)
  
  # Perform FFT with proper shifting
  d0 <- fftshift(fft(fftshift(p1/s1))^n1 * fft(fftshift(p2/s2))^n2)
  d0 <- Re(d0) * exp(n1/(n1+n2) * log(s1) + n2/(n1+n2) * log(s2) - log(k))
  
  # Correct for spacing in frequency domain
  dx <- (b - a) / k
  x <- seq(a, b, by = dx)
  
  data.frame(
    x = x,
    d = d0
  )
}

现在，通过上面的代码，您可以使用这个修改后的函数来计算随机变量之和的密度函数，而不会遇到环绕问题。 这是正态分布的一种示例：

x <- rnorm(100000) + rnorm(100000) + rnorm(100000, -5, 5)
hist(x, freq=FALSE, breaks = 100)
dsum <- dsumf2(dnorm, function(x) dnorm(x, -5, 5), n1 = 2, a = min(x) - 10, b = max(x) + 1)
lines(dsum$x, dsum$d, col = 'red')

这将为您提供三个正态变量之和的正确 PDF，而不会出现周期性问题。

Answer 2

@richarddmoney，您在评论中发布的错误的后续。

您遇到的问题是由傅立叶变换密度及其随后的逆变换的操作引起的。行数不匹配的发生是由于傅里叶域中频率的不正确对齐造成的。

因此，要解决这个问题，需要保证傅里叶变换和后续逆变换所使用的点数的一致性。这需要调整代码以正确处理频域的移位和间隔。

我现在可以整理的修改后的函数和更正如下：

library(SynchWave)  # Load the SynchWave package for fftshift operation

dsumf2 <- function(f1, f2, n1=1, n2=1, a, b, k=2^14) {
  # Perform FFT
  x <- seq(a, b, length.out = k)
  p1 <- f1(x)
  p2 <- f2(x)
  s1 <- sum(p1)
  s2 <- sum(p2)
  
  # Perform FFT with proper shifting
  d0 <- fftshift(fft(fftshift(p1/s1))^n1 * fft(fftshift(p2/s2))^n2)
  d0 <- Re(d0) * exp(n1/(n1+n2) * log(s1) + n2/(n1+n2) * log(s2) - log(k))
  
  # Correct for spacing in frequency domain
  dx <- (b - a) / k
  freq <- fftfreq(k, dx)
  
  # Inverse FFT with proper shifting
  d0 <- ifftshift(d0)
  d0 <- fft(d0, inverse = TRUE)
  d0 <- fftshift(d0)
  
  # Correct for spacing in the time domain
  x <- seq(a, b, by = dx)
  
  data.frame(
    x = x,
    d = Re(d0)
  )
}

在上面更新的函数中，它确保时域中的点数与频域中的点数匹配，从而防止结果数据帧中出现任何不匹配。

现在，您可以使用这个校正后的函数来计算随机变量之和的密度函数，而不会遇到以前的错误。如果您仍然遇到问题，请告诉我。