如何从给定的概率分布中采样？

Question

我正在使用随机数在基态一维框中绘制具有相同概率分布的粒子直方图。但是与原始分布相比，这些值在顶部被削减了。

这里是代码

from numpy import*
from matplotlib.pyplot import*

lit=[]
def f(x):
    return 2*(sin(pi*x)**2)
for i in range(0,10000):
    x=random.uniform(0,1)
    y=random.uniform(0,1)
    if y<f(x):
        lit.append(x)

l=linspace(0,1,10000)
hist(lit,bins=100,density=True)
plot(l,f(l))
show()

生成的图形是：

如何改进此代码以匹配原始代码？

Answer 1

您的问题是您的密度函数 f 在 [0,2] 范围内，但您从 [0,1] 绘制 y。

使用y=np.random.uniform(0,2) 将修复它。

但更好的方法是直接将均匀密度重新映射到您想要的函数：

from scipy.interpolate import interp1d

l = np.linspace(0,1,100)
g = interp1d(l-np.sin(2*np.pi*l)/2/np.pi, l) # inverse for integral of f
# (not sure if analytical solution to above exists... so using interpolation instead)
lit = g(np.random.rand(10000)) 

在f 可能没有解析积分的一般情况下，您也可以使用数值积分：

dl = l[1]
g = interp1d(np.cumsum(f(l))*dl, l)