截断分布与未截断分布不匹配

Question

我想根据截断正态分布生成一个随机列表，但它与其父分布不匹配。

这是我尝试过的代码：

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

# important parameters:
m = 9.1093837e-31    # electron mass (kg)
q = 1.6e-19          # electron charge (C)
k = 1.380649e-23     # boltzmann const (J/K)
phi = 4 * 1.6e-19    # tungsten cathode work function (J)
eps = 8.85e-12       # vacuum permittivity
d = 0.3              # tube length (m)
T_cat = 4000         # cathode temperature (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)

我期望在上面的代码中生成的未截断和截断的分布匹配，但是：

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_in0(T_cat, 0, 100000), bins = 150, histtype = 'step',
         density = True, color = 'red', linewidth=1.2)
plt.show()

任何帮助将不胜感激

Answer 1

听起来您希望两条曲线看起来彼此一致：

您在图中看不到协议，因为：

• 您的

  x_axis

  没有向右延伸足够远，因此您评估未截断正态分布 PDF 的最大值小于截断正态分布的左截断点。
• 截断正态分布经过归一化，使得曲线下的积分为 1，因此它的 y 值将远高于同一 x 坐标处的正态分布。

显示协议的代码以及有关更改内容的注释如下：

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

# important parameters:
m = 9.1093837e-31    # electron mass (kg)
q = 1.6e-19          # electron charge (C)
k = 1.380649e-23     # boltzmann const (J/K)
phi = 4 * 1.6e-19    # tungsten cathode work function (J)
eps = 8.85e-12       # vacuum permittivity
d = 0.3              # tube length (m)
T_cat = 4000         # cathode temperature (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)

# Calculate normalizing constant that relates the magnitudes of the two PDFs
v_min = np.sqrt(2 * phi / m)
std_dev = np.sqrt(k*T_cat/m)
normalizing_constant = 1/norm.sf(v_min/std_dev)

# Extend right limit of `x_axis`
x_axis = np.linspace(1e6, 1e7, 10000)

# Multiply the un-truncated normal PDF by the constant
plt.plot(x_axis, normalizing_constant*norm.pdf(x_axis, 0, np.sqrt(k*T_cat/m)))

# You named the function v_in, not v_in0
plt.hist(v_in(T_cat, 0, 100000), bins = 150, histtype = 'step',
         density = True, color = 'red', linewidth=1.2)

# Zoom in on the part you're interested in
plt.xlim(1.1e6, 1.5e6)
plt.ylim(0, 3e-5)

plt.show()