求助：独立离散随机变量和的分位数高效计算方法

问题背景 #

给定10个独立离散随机变量$X_i$，每个变量有5~20个取值，需计算它们的和$X = X_1 + ... + X_{10}$的分布分位数。直接枚举所有可能的取值组合（总数为$n_1×n_2×...×n_{10}$）规模过大，无法实现。

已尝试两种方法但均遇瓶颈：

• 卷积法：将每个$X_j$近似为$Y_j = X_j + Z$（$Z$是方差极小的正态变量$N(0,σ)$），$Y_j$是多个正态分布的混合，密度振荡极强。$Y=\sum Y_j$的密度是各$Y_j$密度的卷积，但因密度振荡，卷积积分收敛极慢，计算速度无法接受。
• 特征函数法：将$X$近似为$Y=X+Z$，$Y$的特征函数（连续傅里叶变换）可解析计算，但通过逆连续傅里叶变换求解$Y$在数值向量$s$上的密度$f_Y(s)$速度极慢。尝试用逆离散傅里叶变换（IDFT）处理特征函数向量未得到合理结果，且DFT不具备独立变量和的乘积性质，无法直接利用。

疑问 #

逆离散傅里叶变换（IDFT）是否本应能近似逆连续傅里叶变换（ICFT）？

疑问解答 #

IDFT确实可以近似ICFT，但你的尝试失败大概率是没满足以下关键条件：

1. 采样参数匹配：需要选择足够小的频率采样间隔$\Delta\omega$，同时确保采样覆盖的频率范围足够宽，避免截断误差。
2. 周期性修正：IDFT默认处理周期信号，而连续特征函数是非周期的，需用窗函数（如汉宁窗）抑制频谱泄露，或在合适的频率范围截断特征函数，确保误差可控。
3. 尺度因子对应：连续傅里叶变换和离散傅里叶变换的尺度因子必须匹配，否则会出现密度缩放错误。

另外，DFT本身不具备乘积性质，但特征函数的乘积性质依然成立：你可以先计算每个$X_i$的特征函数，相乘得到$X$的特征函数后，再对这个乘积结果做IDFT，而非直接对单个特征函数的采样做DFT再相乘。

可行解决方案 #

针对你的问题，推荐以下高效方法：

1. 改进的特征函数+IDFT方法 #

修正IDFT的使用逻辑，步骤如下：

• 计算每个离散变量$X_i$的特征函数：$cf_{X_i}(t) = \sum_{x} P(X_i=x) e^{itx}$。
• 直接计算和$X$的特征函数：$cf_X(t) = \prod_{i=1}^{10} cf_{X_i}(t)$（无需添加正态噪声$Z$，除非必须平滑分布）。
• 设定频率范围$[-T, T]$与采样点数$N$（建议取2的幂次，加速FFT计算），采样间隔$\Delta t = 2T/N$。
• 对$cf_X(t)$在采样点取值后做IDFT，得到$X$的概率质量函数（PMF）近似（若采样范围覆盖所有可能的$X$取值，可得到精确PMF）。
• 从PMF直接计算分位数。

2. 动态规划剪枝法 #

直接计算$X$的PMF，但通过剪枝忽略概率极小的取值：

• 初始化：以第一个变量$X_1$的PMF作为初始分布。
• 迭代：对每个后续变量$X_i$，将当前PMF与$X_i$的PMF做卷积，丢弃概率小于阈值（如$1e-8$）的取值，大幅减少计算量。
• 最终得到的PMF精度足够支撑分位数计算，且避免了枚举所有组合。

3. 鞍点近似法 #

针对独立变量和的分布，鞍点近似可快速计算分位数，无需精确求解整个分布：

• 计算$X$的累积生成函数（CGF）：$K(t) = \sum_{i=1}^{10} \log\left(\sum_x P(X_i=x) e^{tx}\right)$。
• 求解鞍点方程$K'(t) = q$（$q$为目标分位数对应的取值），得到鞍点$t_0$。
• 利用鞍点近似公式计算累积分布函数（CDF）在$q$处的值，通过二分法找到对应分位数。该方法速度极快，适合大规模变量和的场景。

4. 蒙特卡洛+重要性采样 #

若对精度要求不极致，可采用蒙特卡洛方法结合重要性采样加速：

• 生成大量$X$的样本（每个样本为10个变量取值的和）。
• 直接从样本估计分位数；若样本生成效率低，可通过重要性采样调整权重，重点采样对分位数影响大的区域。

原卷积法问题分析 #

你提供的R代码中，卷积法失效的核心原因是添加的正态噪声$\sigma$过小，导致混合密度振荡过于剧烈，每次卷积都会放大振荡，既慢又容易引入数值误差。若一定要用卷积法，建议：

• 增大$\sigma$的值，平衡平滑效果与精度；
• 使用FFT加速卷积计算，替代数值积分式的卷积。

附原问题中的R代码 #

library(bayesmeta)

# density of (1/10)*sum_{j=1}^10 N(j,0.01) 
# (convex sum of normal distributions)
#
f <- Vectorize(function(s) sum(vapply(1:10,
   FUN = function(j) dnorm(s,mean=j,sd=0.01)/10, FUN.VALUE=0
)))
g <- function(s) dnorm(s,mean=0,sd=0.01)

cat("\n\n")
for(i in 1:5){

  cat("Doing convolution ",i,"\n")
  g <- convolve(g,f)$density
}
cat("\nConvolutions finished, plotting density.")
s <- seq(0,100,length.out=1024)
matplot(s,g(s),type="l")

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