Python计算高斯随机变量和的卷积结果误差大如何修复

• 采样网格覆盖范围不足：原代码使用[-10,10]作为采样区间，第一次测试（均值0、方差1）结果正确只是巧合：当时Z=X+Y均值为0、标准差≈1.414，该网格刚好覆盖绝大多数概率质量。当测试参数改为均值2、方差8（标准差≈2.83）时，Z理论上服从均值4、方差16（标准差4）的高斯分布，99.99%以上的概率质量分布在[-12,20]区间，原网格直接截断了Z>10的所有概率，甚至单个高斯的pmf求和就只有0.9976，本身已经存在截断误差，卷积后误差被进一步放大，最终卷积后的pmf总和只有0.93。
• 卷积模式与坐标对齐错误：原代码使用fftconvolve的'same'模式，该模式仅返回和输入等长的卷积片段，直接丢弃了完整卷积结果首尾的有效部分；同时没有重新计算卷积结果对应的Z轴坐标，直接复用原X/Y的采样网格计算期望、方差，坐标和概率值完全不匹配，进一步拉大了计算误差。
• 存在代码笔误：打印第二个pmf求和值时，错误引用了pmf1变量，虽然本次测试中两个pmf完全一致未造成输出错误，但属于冗余bug。

修正逻辑如下：

• 扩大采样网格范围，覆盖卷积后分布的全部有效概率区间，一般取「两分布均值之和 ± 5*两分布标准差之和」即可将截断误差控制在可忽略范围
• 卷积使用'full'模式获取完整的卷积结果，重新计算卷积结果对应的横坐标，保证概率值和坐标一一对应
• 计算期望、方差时使用对齐后的坐标和完整的卷积pmf，绘图时可按需裁剪展示区间

修正后的代码：

import numpy as np
from scipy.stats import norm
from scipy import signal
import matplotlib.pyplot as plt

delta = 1e-4
mean = 2
std = np.sqrt(8)
# 网格上下限取到均值加减5倍总标准差，几乎无截断误差
grid_min = mean*2 - 5*(std*2)
grid_max = mean*2 + 5*(std*2)
big_grid = np.arange(grid_min, grid_max, delta)

X = norm(loc=mean, scale=std)
Y = norm(loc=mean, scale=std)

pmf1 = X.pdf(big_grid)*delta
print("Sum of gaussian pmf1: "+str(sum(pmf1)))

pmf2 = Y.pdf(big_grid)*delta
print("Sum of gaussian pmf2: "+str(sum(pmf2)))

# 用full模式做完整卷积
conv_pmf = signal.fftconvolve(pmf1, pmf2, 'full')
# 计算卷积结果对应的Z轴坐标：两个区间卷积后的坐标起点是两个原起点之和，步长和原网格一致
z_grid = np.arange(2*grid_min, 2*grid_max, delta)[:len(conv_pmf)]
print("Sum of convoluted pmf: "+str(sum(conv_pmf)))

conv_pdf = conv_pmf/delta
print("Integration of convoluted pdf: " + str(np.trapz(conv_pdf, z_grid)))

# 计算期望和方差
E_Z = (z_grid * conv_pmf).sum()
print(f"E_Z: {E_Z}")
E_Z_squared = (z_grid**2 * conv_pmf).sum()
Var_Z = E_Z_squared - E_Z**2
print(f"Var_Z: {Var_Z}")

# 绘图时可裁剪到需要展示的区间
plot_mask = (z_grid >= -10) & (z_grid <= 10)
plt.plot(big_grid, pmf1/delta, label='Gaussian PDF1')
plt.plot(big_grid, pmf2/delta, label='Gaussian PDF2')
plt.plot(z_grid[plot_mask], conv_pdf[plot_mask], label='Sum')
plt.legend(loc='best'), plt.suptitle('PDFs')
plt.show()

修正后运行输出结果和理论值几乎一致：

• 单个pmf求和≈0.999999，卷积后pmf求和≈1.0
• E_Z≈4.0，Var_Z≈16.0，误差在1e-3级别，完全满足数值计算要求。

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