np.random.dirichlet 小参数：在当前 numpy 中嵌入未来的解决方案

Question

关于当前 np.random.dirichlet 功能的讨论正在进行中，因为它不适用于小参数：

In [1]: import numpy as np

In [2]: np.random.dirichlet(np.ones(3)*.00001)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-2-464b0fe9c6c4> in <module>()
----> 1 np.random.dirichlet(np.ones(3)*.00001)

mtrand.pyx in mtrand.RandomState.dirichlet (numpy/random/mtrand/mtrand.c:25213)()

mtrand.pyx in mtrand.RandomState.dirichlet (numpy/random/mtrand/mtrand.c:25123)()

ZeroDivisionError: float division

可以阅读讨论 here and here 并指出这是规范化错误。目前，出于多种原因，无法将针对小参数切换采样器的建议增强功能合并到 numpy 的主控中。

问题：有人可以建议在 python 中绘制 dirichlets 的不同方法，或者向我指出无需重新编译我的 numpy [=22] 即可使用新采样器的解决方案=] 在未发布的分支上工作？

Answer 1

好的，让我们尝试以下操作。这是 Beta(alpha,beta) 变量采样，它应该适用于任何小数字。

import math
import random

def sample_beta(alpha, beta):
    x = math.log( random.random() )
    y = math.log( random.random() )

    return x / (x + y*alpha/beta)

# some testing
import matplotlib.pyplot as plt

bins = [0.01 * i for i in range(102)]
plt.hist([sample_beta(0.00001, 0.1) for k in range(10000000)], bins)
plt.show()

使用它，您可以尝试按照维基百科中所述通过 Beta 变量对 Dirichlet 进行采样

https://en.wikipedia.org/wiki/Dirichlet_distribution#Random_number_generation

params = [a1, a2, ..., ak]
xs = [sample_beta(params[0], sum(params[1:]))]
for j in range(1,len(params)-1):
    phi = sample_beta(params[j], sum(params[j+1:]))
    xs.append((1-sum(xs)) * phi)
xs.append(1-sum(xs))

如果可行，可以对其进行优化以预先计算所有部分和。

更新

上面的采样依赖于狄利克雷可以通过 beta 变量采样的事实，如果参数较小，这是更好（但更慢）的选择。反过来，beta 变量可以作为一对 gamma 变量进行采样：

beta(a, b) = gamma(1, a) / (gamma(1, a) + gamma(1, b))

如此小的参数从 gamma 中的第一位（如果您直接通过 gamma 变量对 Dirichlet 进行采样）移至第二位。 1（一）在伽马变量中排在第一位意味着它们只是指数分布，采样为 -log(U(0,1))。请检查我的数学是否正确，但这样采样可能有效

np.random.dirichlet 小参数：在当前 numpy 中嵌入未来的解决方案

np.random.dirichlet with small parameter: embed future solution in current numpy

python

random

numpy

dirichlet