时间戳

首页 > 解决方案 > 在 python 中计算自定义概率分布（数字）

python-3.x - 在 python 中计算自定义概率分布（数字）

问题描述

我有一个自定义（离散）概率分布，其形式有点定义：f(x)/(sum(f(x')) for x' in a given离散集X)。此外，0<=x<=1。所以我一直在尝试在 python 3.8.2 中实现它，问题是分子和分母都非常小，python 的浮点表示只是将它们作为 0.0。
在计算完这些概率之后，我需要从一个数组中抽取一个随机元素，它的每个索引可以用分布中的相应概率来选择。所以如果我的分布是 [p1,p2,p3,p4]，我的数组是 [a1,a2,a3,a4]，那么选择 a2 的概率就是 p2，以此类推。
那么我怎样才能以一种优雅而有效的方式实现它呢？
在这种情况下，有什么办法可以使用 np.random.beta() 吗？由于 beta 分布和我的实际分布之间的差异只是归一化常数不同，并且域被限制在几个点。

注意：上面定义的概率质量函数实际上是贝叶斯定理和 f(x)=x^s*(1-x)^f 给出的形式，其中 s 和 f 是给定迭代的固定数。所以确切的问题是，当 s 或 f 变得非常大时，这个东西会变为 0。

标签： python-3.xfloating-pointprecisionbayesianprobability-distribution

解决方案

您可以通过使用日志来计算事物。关键是，虽然分子和分母都可能下溢至 0，但除非您的数字非常小，否则它们的对数不会。

你说

f(x) = x^s*(1-x)^t

所以

logf (x) = s*log(x) + t*log(1-x)

你想计算，比如说

p = f(x) / Sum{ y in X | f(y)}

所以

p = exp( logf(x) - log sum { y in X | f(y)}
  = exp( logf(x) - log sum { y in X | exp( logf( y))}

唯一的困难是计算第二项，但这是一个常见的问题，例如这里

另一方面，手动计算 logsumexp 很容易。

我们想要

S = log( sum{ i | exp(l[i])})

如果 L 是 l[i] 的最大值，则

S = log( exp(L)*sum{ i | exp(l[i]-L)})
  = L + log( sum{ i | exp( l[i]-L)})

最后的总和可以按书面形式计算，因为现在每个项都在 0 和 1 之间，所以没有溢出的危险，并且其中一项（l[i]==L 的项）是 1，所以如果其他术语下溢，即无害。

然而，这可能会失去一点准确性。一种改进是识别索引集 A，其中

l[i]>=L-eps (eps a user set parameter, eg 1)

然后计算

N = Sum{ i in A | exp(l[i]-L)}
B = log1p( Sum{ i not in A | exp(l[i]-L)}/N)
S = L + log( N) + B

推荐阅读