如何生成伪随机卷积？

Question

为了生成伪随机排列，可以使用Knuth 洗牌。对合是一种自逆排列，我想，我可以通过禁止多次触摸一个元素来调整洗牌。然而，我不确定我是否可以有效地做到这一点，以及它是否平等地生成每个对合。

恐怕需要一个例子：在一个集合 上{0,1,2}，有 6 个排列，其中 4 个是对合。我正在寻找一种以相同概率随机生成其中一个的算法。

一个正确但效率很低的算法是：使用 Knuth shuffle，如果没有对合则重试。

Answer 1

我们在这里将其用作a(n)一组大小的对合数n（如 OEIS 所做的那样）。对于给定大小的集合n和该集合中的给定元素，该集合上的对合总数为a(n)。该元素必须在对合过程中保持不变，或者与另一个元素交换。使我们的元素保持固定的对合数量为a(n-1)，因为这些是其他元素的对合。a(n-1)/a(n)因此，对合上的均匀分布必须有保持该元素固定的概率。如果要修复它，只需保留该元素即可。否则，选择另一个尚未被我们的算法检查的元素来与我们的元素交换。我们刚刚决定了集合中的一个或两个元素会发生什么：继续并决定一次会发生一个或两个元素的情况。

为此，我们需要每个 的对合计数列表i <= n，但这可以通过递归公式轻松完成

a(i) = a(i-1) + (i-1) * a(i-2)

（请注意，OEIS 中的这个公式也来自我的算法：第一项计算将第一个元素保持在原处的对合，第二项计算与其交换的元素。）如果您正在使用对合，则此可能足够重要，可以分解为另一个函数，预先计算一些较小的值，并缓存函数的结果以提高速度，如以下代码所示：

# Counts of involutions (self-inverse permutations) for each size
_invo_cnts = [1, 1, 2, 4, 10, 26, 76, 232, 764, 2620, 9496, 35696, 140152]

def invo_count(n):
    """Return the number of involutions of size n and cache the result."""
    for i in range(len(_invo_cnts), n+1):
        _invo_cnts.append(_invo_cnts[i-1] + (i-1) * _invo_cnts[i-2])
    return _invo_cnts[n]

我们还需要一种方法来跟踪尚未决定的元素，以便我们可以有效地以均匀概率选择其中一个元素和/或将元素标记为已决定。我们可以将它们保留在一个缩小的列表中，并在列表的当前末尾添加一个标记。当我们决定一个元素时，我们将当前元素移动到列表末尾以替换已决定的元素，然后减少列表。在这种效率下，该算法的复杂度为O(n)，对每个元素进行一次随机数计算（可能除了最后一个元素）。不可能有更好的订单复杂性。

这是 Python 3.5.2 中的代码。由于通过未定元素列表涉及间接，代码有些复杂。

from random import randrange

def randinvolution(n):
    """Return a random (uniform) involution of size n."""

    # Set up main variables:
    # -- the result so far as a list
    involution = list(range(n))
    # -- the list of indices of unseen (not yet decided) elements.
    #    unseen[0:cntunseen] are unseen/undecided elements, in any order.
    unseen = list(range(n))
    cntunseen = n

    # Make an involution, progressing one or two elements at a time
    while cntunseen > 1:  # if only one element remains, it must be fixed
        # Decide whether current element (index cntunseen-1) is fixed
        if randrange(invo_count(cntunseen)) < invo_count(cntunseen - 1):
            # Leave the current element as fixed and mark it as seen
            cntunseen -= 1
        else:
            # In involution, swap current element with another not yet seen
            idxother = randrange(cntunseen - 1)
            other = unseen[idxother]
            current = unseen[cntunseen - 1]
            involution[current], involution[other] = (
                involution[other], involution[current])
            # Mark both elements as seen by removing from start of unseen[]
            unseen[idxother] = unseen[cntunseen - 2]
            cntunseen -= 2

    return involution

我做了几次测试。这是我用来检查有效性和均匀分布的代码：

def isinvolution(p):
    """Flag if a permutation is an involution."""
    return all(p[p[i]] == i for i in range(len(p)))

# test the validity and uniformness of randinvolution()
n = 4
cnt = 10 ** 6
distr = {}
for j in range(cnt):
    inv = tuple(randinvolution(n))
    assert isinvolution(inv)
    distr[inv] = distr.get(inv, 0) + 1
print('In {} attempts, there were {} random involutions produced,'
    ' with the distribution...'.format(cnt, len(distr)))
for x in sorted(distr):
    print(x, str(distr[x]).rjust(2 + len(str(cnt))))

结果是

In 1000000 attempts, there were 10 random involutions produced, with the distribution...
(0, 1, 2, 3)     99874
(0, 1, 3, 2)    100239
(0, 2, 1, 3)    100118
(0, 3, 2, 1)     99192
(1, 0, 2, 3)     99919
(1, 0, 3, 2)    100304
(2, 1, 0, 3)    100098
(2, 3, 0, 1)    100211
(3, 1, 2, 0)    100091
(3, 2, 1, 0)     99954

对我来说，这看起来非常统一，我检查的其他结果也是如此。