我在一篇论文中编写了Dirichlet-Multinomial后验的推导,我读过这篇文章并且我无法将分布总和加到1.以下是未简化形式的代码:
def pcn(X, n, N, c, alpha):
pnc = np.math.factorial(np.sum([n[i] for i in range(len(n))]))/ \
np.product([np.math.factorial(n[i]) for i in range(len(n))])* \
np.product([c[i]**n[i] for i in range(len(n))])
pc = G(len(X)*alpha)/ \
np.product([G(alpha) for i in range(len(n)) if i in X])* \
np.product([(c[i])**(alpha - 1) for i in range(len(n)) if i in X])
pn = np.math.factorial(N)/ \
np.product([np.math.factorial(n[i]) for i in range(len(n)) if i in X])* \
G(len(X)*alpha)/ \
G(len(X)*alpha + N)* \
np.product([G(alpha + n[i])/G(alpha) for i in range(len(n)) if i in X])
return pnc
我在这里简化了删除分开的部分:
def pcns(X, n, N, c, alpha):
pnc = np.product([c[i]**n[i] for i in range(len(n))])
pc = np.product([(c[i])**(alpha - 1) for i in range(len(n))])/ \
np.product([G(alpha) for i in range(len(n))])
pn = np.product([G(alpha + n[i])/G(alpha) for i in range(len(n))])/ \
G(len(X)*alpha + N)
return pnc * pc / pn
我设置输入变量并初始化c的数组以输入:
X = [0,1]
n = [6.0, 3.0]
N = sum(n)
alpha = 20
c = np.linspace(0, 1, 1000)
然后我迭代c并在每个c和情节评估我的后验函数:
dist = []
for i in c:
dist.append(pcns(X, n, N, [i, 1-i], alpha))
plt.plot(c,dist)
当我总结为dist时得到的值是999或len(c) – 1.有没有人碰巧知道为什么它不是1?
最佳答案 tl; dr:你正在计算一个定积分的离散近似,忘了将这些碎片乘以dx~ = delta_x.
即使这是一个正确归一化的概率分布,为什么dist总和为1?由于我们没有G,我无法准确再现您的结果,所以让我们考虑一个简单的高斯.我们可以将均值设置为0.5,将stddev设置为较小的值,因此0到1的范围应包含大部分概率:
>>> import scipy.stats
>>> c = np.linspace(0, 1, 1000)
>>> p = scipy.stats.norm(0.5,0.02).pdf(c)
>>> sum(p)
999.00000000000045
那又是999.但同样,为什么这应该是1?应该是1的数量是总的综合概率.这里我们只是在某些点上获取概率分布函数的值并将它们加在一起.
更简单的例子:我们知道x ^ 2从0到1的定积分是1/3,但是
>>> x = np.linspace(0, 1, 1000)
>>> sum(x**2)
333.50016683350009
当我们真正想要更像的东西时(粗略的矩形近似,略高于期望的答案,因为我们包括x = 1点的贡献,并且该框实际上超出了积分区域):
>>> sum(x**2) * (x[1]-x[0])
0.3338340008343344
或者在你以前的情况
>>> sum(p) * (c[1]-c[0])
1.0000000000000004
注意,这里我们可以写*(c [1] -c [0]),因为c是等间距的,因此Int中的每个“dx”(p(x)dx)〜= sum(p [x] * delta_x)是一样的.一般来说,我们想要的更像sum(p [: – 1] * np.diff(c)).