截断正态分布样本与非截断正态PDF不匹配的原因及解决方法

问题场景 #

你用以下代码生成截断正态分布的随机变量：

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

# 重要参数：
m = 9.1093837e-31    # 电子质量（kg）
q = 1.6e-19          # 电子电荷量（C）
k = 1.380649e-23     # 玻尔兹曼常数（J/K）
phi = 4 * 1.6e-19    # 钨阴极功函数（J）
eps = 8.85e-12       # 真空介电常数
d = 0.3              # 管长（m）
T_cat = 4000         # 阴极温度（K）

def v_in(T_cat, U, N):
    v_min = np.sqrt(2 * phi / m)
    std_dev = np.sqrt(k*T_cat/m)
    return truncnorm.rvs(v_min/std_dev, np.inf, loc = 0, scale = std_dev, size = N)

预期抽样的直方图能和非截断正态分布的PDF匹配，但运行以下绘图代码后发现二者差异明显：

x_axis = np.linspace(-1e6, 1e6, 10000)
plt.plot(x_axis, norm.pdf(x_axis, 0, np.sqrt(k*T_cat/m)))
plt.hist(v_in(T_cat, 0, 100000), bins = 150, histtype = 'step',
         density = True, color = 'red', linewidth=1.2)
plt.show()

原因分析 #

1. 截断分布的PDF未做归一化：你生成的是原正态分布在v_min右侧的截断分布，该分布的PDF是原正态PDF除以原分布在v_min以上的概率（即生存函数norm.sf(v_min)），而你直接绘制了原正态的PDF，没有做归一化，导致两者高度完全不匹配。
2. x轴范围设置错误：计算可得v_min≈1.185e6 m/s，但你设置的x轴上限只有1e6，既无法展示截断后分布的主要区域，又包含了截断分布不存在的左侧负值区间，视觉上放大了差异。
3. 对截断分布的认知偏差：截断后的分布是原正态分布的子集重新归一化后的结果，它和原分布的形状在截断区间内比例一致，但整体高度必然不同，不可能和原正态PDF完全匹配。

修正方案 #

1. 绘制正确的截断正态PDF #

将原正态PDF除以v_min处的生存函数，只在x≥v_min的区域绘制：

2. 调整x轴范围 #

聚焦到截断分布的有效区间，去掉无意义的左侧区域。

3. 修正后的完整绘图代码 #

# 计算关键参数
v_min = np.sqrt(2 * phi / m)
std_dev = np.sqrt(k*T_cat/m)
# 计算原分布在v_min以上的概率
survival_prob = norm.sf(v_min, loc=0, scale=std_dev)

# 设置合适的x轴范围
x_trunc = np.linspace(v_min, 2.5e6, 10000)
# 绘制截断正态PDF
plt.plot(x_trunc, norm.pdf(x_trunc, 0, std_dev)/survival_prob, label='截断正态PDF')
# 绘制抽样直方图
plt.hist(v_in(T_cat, 0, 100000), bins=150, histtype='step', 
         density=True, color='red', linewidth=1.2, label='抽样直方图')
# 调整x轴显示范围
plt.xlim(v_min*0.9, 2.5e6)
plt.legend()
plt.show()

4. 验证truncnorm参数正确性 #

你的v_in函数中truncnorm.rvs的参数是正确的：a=(v_min - loc)/scale（这里loc=0，scale=std_dev），确保生成的随机变量都大于等于v_min。

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