4

在numpy页面上，他们给出了示例

s = np.random.dirichlet((10, 5, 3), 20)

这一切都很好，很棒；但是如果你想从一个二维的 alpha 数组中生成随机样本呢？

alphas = np.random.randint(10, size=(20, 3))

如果您尝试
np.random.dirichlet(alphas),
np.random.dirichlet([x for x in alphas])或
np.random.dirichlet((x for x in alphas)),

它会导致 ValueError: object too deep for desired array. 唯一似乎有效的是：

y = np.empty(alphas.shape)
for i in xrange(np.alen(alphas)):
    y[i] = np.random.dirichlet(alphas[i])
    print y

...这对于我的代码结构来说远非理想。为什么会这样，谁能想到一种更“类似numpy”的方式来做到这一点？

提前致谢。

4

2 回答 2

4

np.random.dirichlet编写用于为单个 Dirichlet 分布生成样本。该代码是根据 Gamma 分布实现的，并且该实现可以用作矢量化代码的基础，以从不同的分布中生成样本。在下面，dirichlet_sample采用形状为 (n, k) 的数组alphas，其中每一行是alphaDirichlet 分布的向量。它还返回一个形状为 (n, k) 的数组，每一行都是来自 的相应分布的样本alphas。当作为脚本运行时，它会使用dirichlet_sample和np.random.dirichlet来生成样本，以验证它们是否生成相同的样本（直到正常的浮点差异）。

import numpy as np


def dirichlet_sample(alphas):
    """
    Generate samples from an array of alpha distributions.
    """
    r = np.random.standard_gamma(alphas)
    return r / r.sum(-1, keepdims=True)


if __name__ == "__main__":
    alphas = 2 ** np.random.randint(0, 4, size=(6, 3))

    np.random.seed(1234)
    d1 = dirichlet_sample(alphas)
    print "dirichlet_sample:"
    print d1

    np.random.seed(1234)
    d2 = np.empty(alphas.shape)
    for k in range(len(alphas)):
        d2[k] = np.random.dirichlet(alphas[k])
    print "np.random.dirichlet:"
    print d2

    # Compare d1 and d2:
    err = np.abs(d1 - d2).max()
    print "max difference:", err

样品运行：

dirichlet_sample:
[[ 0.38980834  0.4043844   0.20580726]
 [ 0.14076375  0.26906604  0.59017021]
 [ 0.64223074  0.26099934  0.09676991]
 [ 0.21880145  0.33775249  0.44344606]
 [ 0.39879859  0.40984454  0.19135688]
 [ 0.73976425  0.21467288  0.04556287]]
np.random.dirichlet:
[[ 0.38980834  0.4043844   0.20580726]
 [ 0.14076375  0.26906604  0.59017021]
 [ 0.64223074  0.26099934  0.09676991]
 [ 0.21880145  0.33775249  0.44344606]
 [ 0.39879859  0.40984454  0.19135688]
 [ 0.73976425  0.21467288  0.04556287]]
max difference: 5.55111512313e-17
于 2013-04-10T05:01:56.260 回答
2

我想你正在寻找

y = np.array([np.random.dirichlet(x) for x in alphas])

为了您的列表理解。否则你只是传递一个 python 列表或元组。我想numpy.random.dirichlet不接受您的 alpha 值列表的原因是因为它没有设置为 - 它已经接受了一个数组，k根据文档，它期望有一个维度。

于 2013-04-10T01:57:55.040 回答