计算 python 中的自定义概率分布（数值）

Question

我有一个自定义（离散）概率分布，其定义形式如下：f(x)/(sum(f(x')) for x' in a given discrete set 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分布和我的实际分布的区别只是归一化常数不同，域被限制在几个点上。

注：上面定义的Probability Mass函数实际上是贝叶斯定理和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

计算 python 中的自定义概率分布（数值）

Calculating a custom probability distribution in python (numerically)

floating-point

precision

bayesian

python-3.x

probability-distribution