在python中计算自定义概率分布数值

问题描述

我有一个以某种形式定义的自定义（离散）概率分布：给定离散集合X中x'的f（x）/（sum（f（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。

解决方法

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

另一方面，计算logumexp很容易进行。

我们想要

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