概率向量网格化

Question

我正在尝试获取一个n维概率向量网格 - 向量中的每个条目都在0和1之间，所有条目加起来等于1。我希望有可能的每个向量，其中坐标可以取0到1之间均匀间隔的数字v值之一。

为了说明这一点，以下是一个非常低效的实现，当n = 3，v = 3：

from itertools import product
grid_redundant = product([0, .5, 1], repeat=3)
grid = [point for point in grid_redundant if sum(point)==1]

现在grid包含[(0, 0, 1), (0, 0.5, 0.5), (0, 1, 0), (0.5, 0, 0.5), (0.5, 0.5, 0), (1, 0, 0)]。

这种“实现”对于更高维度和更细粒度的网格来说非常糟糕。有没有一种好的方法来解决这个问题，也许使用numpy？

---

我可以添加激励点：如果只是从随机分布中采样给了我足够极端的点，那我可能会很高兴，但事实并非如此。请参见此问题。我所追求的“网格”不是随机的，而是系统地扫描单纯形（概率向量的空间）。

- Schiphol

---

可能是生成总和为1的随机数列表的重复问题。 - yeputons

@yeputons，感谢您的指针。这不是重复的问题；我已经编辑了问题以澄清这一点。 - Schiphol

关于预定概率值，你还能说些什么呢？它们只是[0, 1/(v-1), 2/(v-1), ..., (v-1)/(v-1)]吗？ - Warren Weckesser

是的，抱歉，我想到了均匀间隔的值。我已经编辑了问题以反映这一点。 - Schiphol

Answer 1

以下是递归解决方案。它不使用NumPy，也不是非常高效，但它应该比已发布的代码段更快：

import math
from itertools import permutations

def probability_grid(values, n):
    values = set(values)
    # Check if we can extend the probability distribution with zeros
    with_zero = 0. in values
    values.discard(0.)
    if not values:
        raise StopIteration
    values = list(values)
    for p in _probability_grid_rec(values, n, [], 0.):
        if with_zero:
            # Add necessary zeros
            p += (0.,) * (n - len(p))
        if len(p) == n:
            yield from set(permutations(p))  # faster: more_itertools.distinct_permutations(p)

def _probability_grid_rec(values, n, current, current_sum, eps=1e-10):
    if not values or n <= 0:
        if abs(current_sum - 1.) <= eps:
            yield tuple(current)
    else:
        value, *values = values
        inv = 1. / value
        # Skip this value
        yield from _probability_grid_rec(
            values, n, current, current_sum, eps)
        # Add copies of this value
        precision = round(-math.log10(eps))
        adds = int(round((1. - current_sum) / value, precision))
        for i in range(adds):
            current.append(value)
            current_sum += value
            n -= 1
            yield from _probability_grid_rec(
                values, n, current, current_sum, eps)
        # Remove copies of this value
        if adds > 0:
            del current[-adds:]

print(list(probability_grid([0, 0.5, 1.], 3)))

输出：

[(1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0), (0.5, 0.5, 0.0), (0.0, 0.5, 0.5), (0.5, 0.0, 0.5)]

与已发布的方法快速比较：

from itertools import product

def probability_grid_basic(values, n):
    grid_redundant = product(values, repeat=n)
    return [point for point in grid_redundant if sum(point)==1]

values = [0, 0.25, 1./3., .5, 1]
n = 6
%timeit list(probability_grid(values, n))
1.61 ms ± 20.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit probability_grid_basic(values, n)
6.27 ms ± 186 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

- jdehesa

---

@bobrobbob 嗯，是的，作为一种递归算法，您不需要天文数字大小的输入来打破它...不过我不知道OP期望的实际大小是多少... - jdehesa

抱歉，yup 3000 真的太疯狂了。 - bobrobbob

使用更合理的值，您的结果是错误的。从 n=v=6 开始，您只返回 12 个结果（仅组合为0/.2/1），而“基本”返回 252 个结果（0/.2/.4/.6/.8/1 的组合）。 - bobrobbob

@bobrobbob 谢谢您指出这个问题，我已经修复好了（这是浮点精度误差）。 - jdehesa

Answer 2

在高维向量的完全普遍情况下，即使采用了接受答案中的聪明解决方案，也是相当难以处理的。在我自己的情况下，计算所有值的相关子集是有益的。例如，以下函数计算所有只有n个非零等概率条目的dimension维概率向量：

import itertools as it
import numpy as np

def equip_n(dimension, n):
"""
Calculate all possible <dimension>-dimensional probability vectors with n nonzero,
equiprobable entries
"""
combinations  = np.array([comb for comb in it.combinations(range(dimension), n)])
vectors = np.zeros((combinations.shape[0], dimension))
for line, comb in zip(vectors, combinations):
    line[comb] = 1/n
return vectors 

print(equip_n(6, 3))

这将返回

[[ 0.3333  0.3333  0.3333  0.      0.      0.    ]
 [ 0.3333  0.3333  0.      0.3333  0.      0.    ] 
 [ 0.3333  0.3333  0.      0.      0.3333  0.    ]
 [ 0.3333  0.3333  0.      0.      0.      0.3333]
 [ 0.3333  0.      0.3333  0.3333  0.      0.    ]
 [ 0.3333  0.      0.3333  0.      0.3333  0.    ]
 [ 0.3333  0.      0.3333  0.      0.      0.3333]
 [ 0.3333  0.      0.      0.3333  0.3333  0.    ]
 [ 0.3333  0.      0.      0.3333  0.      0.3333]
 [ 0.3333  0.      0.      0.      0.3333  0.3333]
 [ 0.      0.3333  0.3333  0.3333  0.      0.    ]
 [ 0.      0.3333  0.3333  0.      0.3333  0.    ]
 [ 0.      0.3333  0.3333  0.      0.      0.3333]
 [ 0.      0.3333  0.      0.3333  0.3333  0.    ]
 [ 0.      0.3333  0.      0.3333  0.      0.3333]
 [ 0.      0.3333  0.      0.      0.3333  0.3333]
 [ 0.      0.      0.3333  0.3333  0.3333  0.    ]
 [ 0.      0.      0.3333  0.3333  0.      0.3333]
 [ 0.      0.      0.3333  0.      0.3333  0.3333]
 [ 0.      0.      0.      0.3333  0.3333  0.3333]]

这非常快。 %timeit equip_n(6, 3) 返回

15.1 µs ± 74.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

- Schiphol

---

这不是一个解决方案。首先，如果v = 3，则应该从0、0.5和1中选择概率值。其次，它并不包括所有可能的向量，而只包括其中一些。例如[1, 0, 0, 0, 0, 0]。 - gciriani