Two cycles.py
From Werner KRAUTH
| Revision as of 20:36, 5 March 2025 Werner (Talk | contribs) ← Previous diff |
Revision as of 21:40, 5 March 2025 Werner (Talk | contribs) Next diff → |
||
| Line 3: | Line 3: | ||
| __FORCETOC__ | __FORCETOC__ | ||
| =Description= | =Description= | ||
| + | |||
| + | To sample permutations of N elements which only contain cycles of length 1 or 2, we derived in Lecture 7 a recursion formula for the number Y_N of such permutations (eq. (7.6)), namely | ||
| + | Y_N = Y_{N-1} + (N - 1) Y_{N-2} | ||
Revision as of 21:40, 5 March 2025
This is a Python3 program to sample random permutations P of N elements subject to the constraint that the lengths of the cycles in P can only be one or two. The algorithm is based on a recursive-sampling strategy discussed in my book and presented in Lecture 7 of my 2025 Oxford lectures.
Contents |
Description
To sample permutations of N elements which only contain cycles of length 1 or 2, we derived in Lecture 7 a recursion formula for the number Y_N of such permutations (eq. (7.6)), namely
Y_N = Y_{N-1} + (N - 1) Y_{N-2}
Program
import random
from sympy.combinatorics import Permutation
Y = {-1: 0, 0:1, 1:1}
N = 4
Stats = {}
for k in range(2, N + 1):
Y[k] = Y[k - 1] + (k - 1) * Y[k - 2]
for iter in range(100000):
Q = list(range(N))
random.shuffle(Q)
M = N
P = []
while M > 0:
if random.uniform(0.0, Y[M]) < Y[M - 1]:
P.append([Q[M - 1]])
M -= 1
else:
P.append([Q[M - 1] , Q[M - 2]])
M -= 2
P = tuple(Permutation(P).array_form)
Stats[P] = Stats.get(P, 0) + 1
print(Stats)
Output
The output is a dictionary of permutations (in array form) together with the number of times they were sampled. Here we see, that all the permutations are (roughly) equally frequent, that is, were drawn with equal probabilities. We may check that all these permutations have cycles of at most two.
{(0, 3, 2, 1): 9801, (0, 2, 1, 3): 10281, (0, 1, 2, 3): 10145, (2, 3, 0, 1): 9940, (3, 1, 2, 0): 9869,
(0, 1, 3, 2): 9950, (2, 1, 0, 3): 9830, (3, 2, 1, 0): 10154, (1, 0, 3, 2): 10016, (1, 0, 2, 3): 10014}
Version
See history for version information. For more information on this program, see Lecture 7 of my 2025 Oxford lectures.
