Probability that Secret Santa arrangement will result in perfect pairings

So, we had Secret Santa at work.

We’re 8 people.
We each took turns and pulled a small piece of paper from a bowl with a name on it.
The only rule : If you pull your name, you have to put the piece of paper back in the bowl and try again.

Let’s call the people A, B, C, D, E, F, G, H, which is also the order in which they picked their piece of paper.

We did the gift exchange last night.

A was F’s secret santa.
B was E’s secret santa.
C was D’s secret santa.
D was C’s secret santa.
E was B’s secret santa.
F was A’s secret santa.
G was H’s secret santa.
H was G’s secret santa.

See what happened? We made couples.

A and F were each other’s secret santa.
B and E were each other’s secret santa.
C and D were each other’s secret santa.
G and H were each other’s secret santa.

What is the probability of this happening and how do you calculate it?

The total number of assignments among $2n$ people, where nobody is assigned to themselves, is (These are called derangements.) The value is very close to $(2n)! / e$.

If they correspond to perfect pairings, then they are a product of disjoint transpositions. This implies their cycle structure is of the form

The number of distinct such patterns is the order of the group of all permutations of the $2n$ names divided by the order of the stabilizer of the pattern. A stabilizing element may swap any number of the pairs and it may also permute the $n!$ pairs, whence there are $2^n n!$ stabilizing elements. Therefore there are

such pairings.

Since all such perfect pairings are derangements, and all derangements are equally likely, the chance equals

For $2n=8$ people the exact answer therefore is $15/2119 \approx 0.00707881$ while the approximation is $e/(2^4\, 4!) \approx 0.00707886$: they agree to five significant figures.

To check, this R simulation draws a million random permutations of eight objects, retains only those that are derangements, and counts those that are perfect pairings. It outputs its estimate, the standard error of the estimate, and a Z-score to compare it to the theoretical value. Its output is

       p.hat           se            Z
0.006981031  0.000137385 -0.711721705


The small Z-score is consistent with the theoretical value. (These results would be consistent with any theoretical value between $0.0066$ and $0.0073$.)

paired <- function(x) crossprod(x[x] - 1:length(x))==0
good <- function(x) sum(x==1:length(x)) == 0

n <- 8
set.seed(17)
x <- replicate(1e6, sample(1:n, n))
i.good <- apply(x, 2, good)
i.paired <- apply(x, 2, paired)

n.deranged <- sum(i.good)
k.paired <- sum(i.good & i.paired)
p.hat <- k.paired / n.deranged
se <- sqrt(p.hat * (1-p.hat) / n.deranged)
(c(p.hat=p.hat, se=se, Z=(p.hat - 15/2119)/se))