- android - 多次调用 OnPrimaryClipChangedListener
- android - 无法更新 RecyclerView 中的 TextView 字段
- android.database.CursorIndexOutOfBoundsException : Index 0 requested, 光标大小为 0
- android - 使用 AppCompat 时,我们是否需要明确指定其 UI 组件(Spinner、EditText)颜色
我正在对一篇我读过的论文中 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
sum 应该为 1?由于我们没有 G
我无法准确重现您的结果,所以让我们考虑一个简单的高斯。我们可以将均值设置为 0.5 并将标准差设置为较小的值,这样 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?应为一的数量是总的综合概率。这里我们只是取概率分布函数在某些点的值并将它们相加。
更简单的例子:我们知道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
是等间距的所以中的每个“dx” >Int(p(x) dx) ~= sum(p[x] * delta_x)
是一样的。一般来说,我们想要更像 sum(p[:-1] * np.diff(c))
的东西。
关于python - 在 Python 中规范化后验分布时遇到问题,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/33852172/
我是一名优秀的程序员,十分优秀!