为了生成伪随机排列，可以使用 Knuth shuffles。对合是自逆排列，我想，我可以通过禁止多次触摸元素来适应洗牌。但是，我不确定我是否能有效地做到这一点，以及它是否能等价地产生每一次对合。
恐怕需要一个例子:在集合 {0,1,2} 上，有 6 个排列，其中 4 个是对合。我正在寻找一种算法，以相同的概率随机生成其中一个。
一个正确但非常低效的算法是:使用 Knuth shuffle，如果没有对合则重试。
最佳答案
我们在这里使用 a(n) 作为一组大小为 n ( as OEIS does ) 的对合次数。对于给定大小的集合 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
这对我来说看起来很统一，我检查的其他结果也是如此。
关于algorithm - 如何生成伪随机对合？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/38995387/