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

Question

我有一个自定义的（离散的）概率分布，在某种程度上定义为：f(x)/(sum(f(x')) for x' in a given离散集X)。此外，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。

Answer 1

您可以通过使用日志来很好地计算事物。关键是，虽然分子和分母都可能下溢为 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))}

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

另一方面，手动计算 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