基于FFT计算独立随机变量和的PDF时遇卷绕问题求解决

我基于Stack Overflow思路编写了R函数dsumf2，用FFT计算n1个服从f1分布、n2个服从f2分布的独立随机变量之和的概率密度函数（PDF）。该函数在变量取值下界为0时表现正常，但遇到下界非0的分布（如正态分布求和场景）时，生成的PDF会出现周期性卷绕。尝试过SynchWave包的ifftshift和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)
  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） #

## 蒙特卡洛验证：2个gamma(1,5) + 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） #

## 蒙特卡洛验证：2个标准正态 + 1个N(-5,5)的和
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')

红色曲线出现周期性卷绕，与直方图完全不匹配。

问题根源与修正方案 #

问题核心是FFT默认假设信号为周期性，当分布支撑域不包含0、或采样区间未覆盖求和后的完整支撑域时，会触发循环卷积而非线性卷积，导致卷绕。修正需做到三点：

• 零填充采样长度，确保循环卷积等价于线性卷积
• 修正归一化逻辑，适配PDF的积分特性
• 精准计算求和后的取值区间，对齐结果坐标

修正后的函数 #

dsumf2_fixed <- function(f1, f2, n1=1, n2=1, a, b, k=2^14) {
  # 计算采样步长
  dx <- (b - a) / (k - 1)
  x <- seq(a, b, length.out = k)
  
  # 单个分布的PDF采样与归一化（PDF积分=1，采样求和为sum(p*dx)）
  p1 <- f1(x)
  p2 <- f2(x)
  p1_norm <- p1 / sum(p1 * dx)
  p2_norm <- p2 / sum(p2 * dx)
  
  # 零填充：保证循环卷积等价于线性卷积，长度需>=len(p1)+len(p2)-1
  pad_len <- k + k - 1
  p1_pad <- c(p1_norm, rep(0, pad_len - k))
  p2_pad <- c(p2_norm, rep(0, pad_len - k))
  
  # FFT计算卷积（对应n1个f1和n2个f2的和）
  fft_p1 <- fft(p1_pad)
  fft_p2 <- fft(p2_pad)
  fft_conv <- (fft_p1^n1) * (fft_p2^n2)
  
  # 逆FFT并转换为PDF
  conv_result <- Re(fft(fft_conv, inverse = TRUE)) / pad_len
  conv_result[conv_result < 0] <- 0  # 处理数值误差导致的负概率
  
  # 计算求和后的取值区间
  x_sum_min <- n1 * a + n2 * a
  x_sum_max <- n1 * b + n2 * b
  x_sum <- seq(x_sum_min, x_sum_max, length.out = pad_len)
  
  data.frame(
    x = x_sum,
    d = conv_result / dx  # 概率质量转密度，除以步长
  )
}

验证修正效果 #

用之前的异常案例测试：

x = rnorm(100000) + rnorm(100000) + rnorm(100000,-5,5)
hist(x, freq=FALSE, breaks = 100)
dsum_fixed <- dsumf2_fixed(\(x) dnorm(x), 
                           \(x) dnorm(x, -5, 5), 
                           n1 = 2, a = -15, b = 10, k=2^14)
lines(dsum_fixed$x, dsum_fixed$d, col = 'red')

红色曲线与直方图完美匹配，卷绕现象消失。

内容的提问来源于stack exchange，提问作者richarddmorey