为啥从我的自定义分布中抽取的随机样本不遵循 pdf？

Posted 2023-03-12

【中文标题】为啥从我的自定义分布中抽取的随机样本不遵循 pdf？

【英文标题】：Why are the random samples drawn from my custom distribution not following the pdf?为什么从我的自定义分布中抽取的随机样本不遵循 pdf？

【发布时间】：2020-07-03 01:06:50

【问题描述】：

我使用 scipy 的 rv_continuous 方法创建了一个自定义发行版。我正在尝试创建由 beta 衰变产生的电子的能量分布。鉴于其pdf：

我来自：http://hyperphysics.phy-astr.gsu.edu/hbase/Nuclear/beta2.html#c1

我定义我的分布：

import numpy as np
from scipy.stats import rv_continuous 
import matplotlib.pyplot as plt 

class beta_decay(rv_continuous):
    def _pdf(self, x):
        return (22.48949986*np.sqrt(x**2 + 2*x*0.511)*((0.6-x)**2)*(x+0.511))

# create distribution from 0 --> Q value = 0.6 
beta = beta_decay(a=0, b= 0.6)

# plot pdf 
x = np.linspace(0,0.6)
plt.plot(x, beta.pdf(x))
plt.show()

# random sample the distribution and plot histogram 
random = beta.rvs(size =100)
plt.hist(random)
plt.show()

其中 x = KE, Q = 0.6, C = 22.48...（通过在 0 --> Q 和设置等于 1 之间整合上述表达式来归一化），并且我忽略了费米函数 F(Z' ,KEe) 在上式中。

当我绘制 pdf 时，它看起来正确：

但是，当我尝试使用 .rvs() 从中抽取随机样本时，它们所取的值在 RHS 处大幅达到峰值，而不是像我预期的那样低于 pdf 的峰值：

最终，我的代码需要对分布进行采样，以获得由 beta 衰变释放的电子的 KE。为什么我的直方图如此错误？

【问题讨论】：

请将代码发布为文本，而不是屏幕截图。

在我看来您的 PDF 是错误的。仔细观察，该图中曲线下的面积看起来比 1 小很多，并且您的 _pdf 方法中的表达式将迅速增长到 x=0.6 以上。

"并且我忽略了上述方程中的费米函数 F(Z',KEe)" - 你为什么这样做？

@user2357112supportsMonica 您好，现在以代码形式发布 - 抱歉！我将尝试添加费米函数——我的理解是它只是一个校正因子——在这一点上，我只想要一些有效的东西（即我可以正确定义分布并从中采样）而不是物理精度。我只希望我的 pdf 定义在 0 -> 0.6 之间，因为 0.6 是放射性衰变可用的最大能量。也许我对 pdf 的规范化不正确。

【参考方案1】：

我认为您的 PDF 以错误的方式定义，它没有被规范化。在我对其进行归一化并制作正确的直方图后，它似乎工作正常

代码（Win10 x64，Anaconda Python 3.7）

#%%
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

from scipy.stats import rv_continuous

def bd(x):
    return (22.48949986*np.sqrt(x**2 + 2*x*0.511)*((0.6-x)**2)*(x+0.511))

a = 0.0
b = 0.6

norm = integrate.quad(bd, a, b) # normalization integral
print(norm)

class beta_decay(rv_continuous):
    def _pdf(self, x):
        return bd(x)/norm[0]

# create Q distribution in the [0...0.6] interval
beta = beta_decay(a = a, b = b)

# plot pdf
x = np.linspace(a, b)
plt.plot(x, beta.pdf(x))
plt.show()

# sample from pdf
r = beta.rvs(size = 10000)
plt.hist(r, range=(a, b), density=True)
plt.show()

还有情节

抽样

【讨论】：

谢谢！这正是我一直在寻找的。而且您向我展示了与我不知道存在的 scipy 集成的整个世界：D。接受你的回答。