利用傅里叶变换验证指数和服从Gamma分布的代码问题排查

我想利用结论——Exp(β)分布的傅里叶变换$p(\omega)$的n次方对应Gamma(n, β)分布——编写代码验证指数和的分布特性，已设β=1简化计算。代码在n较小时效果正常，但n增大后Gamma分布曲线在定义域上偏移明显，找不到代码错误。注：掩码用于消除数值傅里叶变换后常出现的水平线伪影。

以下是代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import expon, gamma

# 确定性场景
n = 15   # 傅里叶变换的缩放因子
N = 1024     # 采样点数
T = 0.01     # 采样间隔
shift_point = 100 # 期望的图像中心位置

# 离散化PDF
xm = 0  # 指数分布的定义域起点
x_values = np.linspace(xm, 100, N, endpoint=False)
pdf_values = expon.pdf(x_values, scale=1)

# 指数分布PDF的数值傅里叶变换
P_omega = np.fft.fft(pdf_values)
# 计算n次幂，对应n个独立指数变量和的特征函数
exp_nP = P_omega**(n)

# 构造用于逆傅里叶变换的完整频率数组
full_exp_nP = np.concatenate([exp_nP, np.conj(exp_nP[::-1])])

# 逆傅里叶变换
ift_result = np.fft.ifft(full_exp_nP)
ift_result = ift_result[::2]

# 绘图预处理：移除小值消除伪影
mask = ift_result > 0.01
ift_result = ift_result[mask]
x = np.fft.fftfreq(N, T)[mask]

plt.figure(figsize=(10, 5))
plt.plot(x, -np.abs(ift_result)/np.trapz(x, np.abs(ift_result)), label = "傅里叶变换后逆变换结果")
plt.title('Exp(n·P(ω))的逆傅里叶变换（中心对齐100）')
plt.xlabel('x')
plt.ylabel('概率密度')
plt.grid(True)

# 绘制Gamma分布作为理论对比
x_gamma = np.arange(0, 50, 0.01)
y_gamma = gamma.pdf(x_gamma, n, scale=1)
plt.plot(x_gamma, y_gamma, label = "Gamma(n, 1)分布")
plt.legend()
plt.show()

图像对比 #

• n=10时：傅里叶逆变换结果与Gamma曲线吻合较好
• n=30时：傅里叶逆变换结果与Gamma曲线偏差明显，出现定义域偏移

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