0

关于当前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

可以在此处和此处阅读讨论并指出这是一个规范化错误。目前，出于多种原因，针对小参数切换采样器的提议增强无法合并到 numpy 的 master 中。

问题：有人可以建议一种不同的方法来在 python 中绘制 dirichlets 或者指出一个解决方案来使用新的采样器而不重新编译我的 numpy 和/或在未发布的分支上工作？

4

1 回答 1

2

好的，让我们尝试以下方法。这是 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))

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

更新

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

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

因此，小参数从 gamma 中的第一个（如果您直接通过 gamma 变量对 Dirichlet 进行采样）移动到第二个。并且 1（一）在 gamma 变量中排在第一位意味着它们只是指数分布，采样为 -log(U(0,1))。请检查我的数学是否正常，但这样采样可能会起作用

于 2015-11-13T01:38:24.177 回答