基于FFT求解高斯分布之和PDF时的奇偶数量结果差异问题

嘿，这个问题挺有意思的——按说不管加多少个独立高斯随机变量，结果都得是高斯分布啊，毕竟高斯分布的可加性是它的核心性质之一。你碰到的这种奇偶数量结果不一样的情况，大概率是代码里某个细节跑偏了，咱们一起来排查排查：

可能的问题原因 #

1. 不小心引入了绝对值/截断操作 #

如果你的代码里，在计算单个变量的PDF、或者做卷积的过程中，对偶数次求和的结果误加了绝对值，或者写了类似“当变量数为偶数时只保留非负区间”的条件判断，那肯定会得到两个半高斯拼接的结果。比如要是有if n % 2 == 0: pdf = np.where(x >= 0, pdf, 0)这种代码，直接就把负半轴的分布砍掉了，看起来自然不对。

2. 卷积实现的逻辑bug #

很多人手动实现卷积的时候容易在边界处理上踩坑：

• 比如用np.convolve时选了mode='same'，这个模式会截断卷积后的边缘结果，次数多了（尤其是偶数次）就会把负半轴的分布信息弄丢，最后看起来像半高斯拼接。正确的做法应该用mode='full'来保留完整的卷积结果范围。
• 或者循环卷积的时候，每次没有正确更新x轴的取值范围，导致多次卷积后x轴范围不够，分布的尾部被截断，也会出现类似的奇怪形状。

3. 采样/数值计算的偏差 #

比如生成单个高斯变量的时候，是不是误生成了只有非负的样本？或者计算PDF时设置的x轴网格范围太小，偶数次求和后分布的方差是n*σ²，比奇数的更大，结果边缘被截断，看起来就像两个半高斯拼起来了。

给你一个正确的参考实现 #

你可以对照这个代码看看自己的问题出在哪：

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def sum_gaussian_pdf(n, mean=0, std=1):
    # 理论上，n个独立高斯之和的分布是N(n*mean, sqrt(n)*std)
    # 先设置足够宽的x轴范围，避免截断
    x = np.linspace(-5 * np.sqrt(n), 5 * np.sqrt(n), 1000)
    theoretical_pdf = stats.norm.pdf(x, loc=n*mean, scale=std*np.sqrt(n))
    
    # 用数值卷积模拟求和过程
    single_pdf = stats.norm.pdf(x, loc=mean, scale=std)
    current_pdf = single_pdf.copy()
    for _ in range(n-1):
        # 用full模式做卷积，保留完整结果
        current_pdf = np.convolve(current_pdf, single_pdf, mode='full')
        # 更新x轴范围：两次卷积后x的范围是原来的两倍
        new_x = np.linspace(x[0]*2, x[-1]*2, len(current_pdf))
        x = new_x
        # 归一化，保证PDF积分是1
        current_pdf /= np.trapz(current_pdf, x)
    
    return x, current_pdf, theoretical_pdf

# 测试偶数和奇数情况
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# 偶数n=4
x4, pdf4, theo4 = sum_gaussian_pdf(4)
ax1.plot(x4, pdf4, label='数值卷积结果')
ax1.plot(x4, theo4, label='理论PDF', linestyle='--', color='orange')
ax1.set_title('4个独立高斯变量之和')
ax1.legend()

# 奇数n=3
x3, pdf3, theo3 = sum_gaussian_pdf(3)
ax2.plot(x3, pdf3, label='数值卷积结果')
ax2.plot(x3, theo3, label='理论PDF', linestyle='--', color='orange')
ax2.set_title('3个独立高斯变量之和')
ax2.legend()

plt.tight_layout()
plt.show()

这个代码不管n是奇数还是偶数，数值卷积的结果都会和理论高斯分布完全重合。如果你的代码运行结果还是不对，可以把你的代码片段贴出来，咱们再仔细找问题~

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