python - 为什么从我的自定义分布中抽取的随机样本不遵循 pdf?
问题描述
我使用 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)在上面的方程式中。
但是,当我尝试使用 .rvs() 从中抽取随机样本时,它们所取的值会在 RHS 上达到峰值,而不是像我预期的那样低于 pdf 的峰值:
最终,我的代码需要对分布进行采样以获得由 beta 衰变释放的电子的 KE。为什么我的直方图如此错误?
解决方案
我认为您的 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()
和情节
采样
推荐阅读
- java - 如何将两个活动制作到一个文件中?
- f# - f# 将 seq obj 转换为 seq 记录
- javascript - 使用 Content-Type 时 nodejs https 请求抛出错误
- ruby-on-rails - Heroku 上的 WickedPDF 松散样式
- javascript - 是否可以在我可以在同一个 HTML 文件中编辑的 HTML 文件中创建一个 Javascript 对象?
- javascript - 对 XMLHttpRequest 使用“GET”时,如何限制传入信息的大小?
- php - 下载特定文件 URL 时触发操作
- javascript - 从 TypeScript 中的每个组中检索 n 个项目
- sql - 使用内联查询从 SSRS 数据集中提取数据库对象依赖项(查询类型:文本)
- mysql - 如何进行查询以显示单个递归请求的结果?