Probability that Secret Santa arrangement will result in perfect pairings

Question

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?

Answer 1

The total number of assignments among $2n$ people, where nobody is assigned to themselves, is $$d(2n) = (2n)!(1/2 - 1/6 + \cdots + (-1)^k/k! + \cdots + 1/(2n)!).$$ (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

$$(a_{11}a_{12})(a_{21}a_{22})\cdots(a_{n1}a_{n2}).$$

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

$$p(2n) = \frac{(2n)!}{2^n n!}$$

such pairings.

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

$$\frac{p(2n)}{d(2n)} = \frac{1}{2^n n!(1 - 1/2 + 1/6 - \cdots + (-1)^k/k! + \cdots + 1/(2n)!)} \approx \frac{e}{2^n n!}.$$

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))

Answer 2

I was quite impressed by the elegance in @whuber answer. To be honest I had to do a lot of acquainting myself with new concepts to follow the steps in his solution. After spending a lot of time on it, I've decided to post what I got. So what follows is an exegetical note to his already accepted response. In this way there is no attempt at originality, and my only objective is to provide some additional anchoring points to follow some of the steps involved.

So here it goes...

1. Why $2n$? Well this one may be way too basic: we need an even number of people to play the game.

2. Can we derive the formula for derangements?

Following the Wikipedia entry with the example based on matching of people and hats, $n$ hats to be precise, we have that the options for the first person picking up a hat are shown here as follows:

$d(n)=(n-1)[d(n-2) + d(n-1)]=$

$=n\,d(n-2)-d(n-2)+n\,d(n-1)-d(n-1)$, which can be reorganized as:

$d(n) - n\,d(n-1) = -[d(n-1)-(n-1)\,d(n-2)]$.

Now noticing the parallelism between the LHS of this equation and the part on the RHS within brackets we can continue recursively:

$d(n) - n\,d(n-1)=-[d(n-1)-(n-1)\,d(n-2)] =$

$=(-1)^2\,[d(n-2)-(n-2)\,d(n-3)]=\cdots=(-1)^{n-2}\,d(2)-2\,d(1)$

This implies that $d(n) = n\,d(n-1)+(-1)^n$.

Working backwards:

$d(2)=1$

$d(3)= 3\,d(2)-1 =3*1\,-1$

$d(4)= 4\,d(3)+1 =4*3*1\,-4\,\,+1$

$d(5)= 5\,d(4)-1 =5*4*3*1\,-5*4\,+5\,\,-1$

$d(6)= 6\,d(5)+1 = 6*5*4*3*1\,-6*5*4\,+6*5\,-6\,+1=$

$=6!\left(\frac{1}{2}- \frac{1}{3*2}+\frac{1}{4*3*2}-\frac{1}{5*4*3*2}+\frac{1}{6!}\right)=$

$=6!\left(\frac{1}{6!}-\frac{1}{5!}+\frac{1}{4!}-\frac{1}{3!}+\frac{1}{2!}-\frac{1}{1!}+1\right)$

So in general,

$d(n)=n!\left(1-1+\frac{1}{2!}-\frac{1}{3!}+\frac{1}{4!}+\cdots+\frac{1}{n!}\right)$

And remembering the Taylor series of $e^x$ evaluated at $x=-1$:

$d(n)\approx\frac{n!}{e}$

3. Disjoint transposition: The concept of transposition is easy to get from the link provided in the original answer, but "disjoint" was a bit less clear. Looking at an example, from the set ${a,b,c,d,e,f}$, the permutation ${b,d,a,c,f,e}$ can be expressed as a cycle as: a -> b -> d -> c after which it returns to a, but e -> f forms a loop onto itself - a disjoint cycle. Or, putting these two cycles together, the permutation can be expressed as the product $(\text {a b d c})(\text{e f})$.

In the Santa Claus question (eight employees happening to have drawn names in perfectly matched pairs: Anna gives Martha a gift, while Martha has drawn Anna's name) there will be $4$ closed loops.

4. To find the number of two-element loops we need to divide all possible permutations $(2n)!$ of the set of eight ($2n$) people by the total possible number of swaps of two elements $2^n$ and the total number of permutations of these pairs $n!$: $p(2n) = \frac{(2n)!}{2^n n!}$.

---

For the R simulation:

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

This function boils down to understanding x[x]: It is meant to evaluate an $8$ element vector representing the present assignments, and determine whether it is composed of 2-element loops, as in the Santa Claus problem. As long as the permutations correspond to transposing elements so that if Paul was supposed to give a present to Maria (Paul -> Maria and vice versa, Maria -> Paul) and Max to John (Max -> John / John -> Max) initially, the resultant transposition just results in new perfect pairing (Max -> Maria / Maria -> Max and Paul -> John / John -> Paul) we are fulfilling the initial condition of perfect pairing:

In other words, if we go back to the example of the hats in Wikipedia person i always takes back hat $1$.

2. good <- function(x) sum(x==1:length(x)) == 0

This function evaluates whether we are dealing with a derangement by comparing the vector $x$ element wise to the the vector $(1,2,3,4,5,6,7,8)$, and making sure there is no coincidence.

3. k.paired <- sum(i.good & i.paired) is there to exclude paired permutations like the one above in the diagram, which are not derangements:

v <- c(1,2,3,4,5,6,7,8)
w <- c(1,2,3,5,4,6,7,8)

(c("is v paired?" = paired(v), "is w paired?" = paired(w),
   "is v a derang?" = good(w), "is w a derang?" = good(w)))

# not all paired permutations are derangements.

Answer 3

So this answer is prompted by a few questions in the comments. The earlier two answers are correct if the implicit assumption that all derangements are equally likely is true. This is sadly not the case with the method actually used in the question which is why quite a lot of secret santas are now using computer based methods as these can afford to discard the full drawing and re-start when the random assignment gives someone themselves.

To check that the difference actually was real, I implemented a set of python functions to systematically generate all possible options for a given set of people.

Firstly I implemented a version using the method that a computer will use which is generate a completely random option and repeat till we get a valid answer. This gives the same result as just counting the number of valid options (which is the same as the number of derangements) and assigning equal probability to each of them.

Second added in the handling of the replace and re-draw if you select your own name. As written a person could keep drawing them selves, putting their name back in and repeating for an infinite number of times, however we can show that this is actually the same as drawing any of the valid options with equal probability as $$ \frac{1}{n}\sum_{i=0}^\infty\left(\frac{1}{n}\right)^i= \frac{1}{n-1} $$

This gives us another way of viewing the problem. The probability of selecting each option now depends on if your name is still in the hat or not. We aren't suggesting people should check if they are in the hat, just that the probabilities are the same. The effect this difference in probabilities has is also compounded by the effect of the having to still abandon and re-start if we end up with the last name in the hat being the name of the last person to pick.

Sadly I haven't found an explicit formula for the probabilities of each option but I have implemented python code with this calculation and produced graphs of the resulting probabilities and calculated the probabilities for the pairing to happen.

The calculation done in the other answers gives a probability of 0.0070788107597923545 where as the accurate calculation for the stated method gives a probability of 0.006876714828601176. This isn't a large difference but it definitely is a measurable difference.

My code is included below for reference:

import pylab as pl

def generate_options(count):
    entries = set(range(count))
    return do_generate_options([], entries)

def do_generate_options(selected, available):
    """ Calculate the probabilities of a secret santa draw if we abandom if someone draws their own name. """
    if len(available) == 0:
        yield 1,tuple(selected)
        return
    factor = 1.0/len(available)
    for i in list(available):
        if i == len(selected):
            yield factor, None
        else:
            selected.append(i)
            available.remove(i)
            for p, option in do_generate_options(selected, available):
                yield p * factor, option
            selected.pop()
            available.add(i)

def do_generate_options(selected, available):
    """ Calculate the probabilities of a secret santa draw if you return the name and re-draw from that person if someone draws their own name. """
    if len(available) == 0:
        yield 1,tuple(selected)
        return
    factor = 1.0/len(available)
    if len(selected) in available:
        if len(available) == 1:
            yield 1, None
            return
        # sum 1/n^i for i from 1 to infinity = 1/n + 1/n * 1/(n-1) = 1/(n-1)
        factor = 1.0/(len(available) - 1)
    for i in list(available):
        if i != len(selected):
            selected.append(i)
            available.remove(i)
            for p, option in do_generate_options(selected, available):
                yield p * factor, option
            selected.pop()
            available.add(i)

def plot_probability_counts(count):
  data = {}
  for p, entry in generate_options(count):
    if entry is not None:
      if p in data:
        data[p] = data[p] + 1
      else:
        data[p] = 1
  entries = list(data.items())
  plt = pl.figure().add_subplot(111)
  plt.set_xlabel("probability")
  plt.set_ylabel("count")
  plt.plot([p * total for p, c in entries], [c for p, c in entries], "+")

def plot_probability_counts_scaled(count):
  data = {}
  total = 0
  for p, entry in generate_options(count):
    if entry is not None:
      total = total + 1
      if p in data:
        data[p] = data[p] + 1
      else:
        data[p] = 1
  entries = list(data.items())
  plt = pl.figure().add_subplot(111)
  plt.set_xlabel("probability/uniform probability")
  plt.set_ylabel("fraction")
  plt.plot([p * total for p, c in entries], [c / total for p, c in entries], "+")

def plot_probabilities_per_person(count):
  total = 0
  probabilites = [[0] * count for i in range(count)]
  for p, entry in generate_options(count):
    if entry is not None:
      total = total + p
      for i,j in enumerate(entry):
        probabilites[i][j] = probabilites[i][j] + p
  width = 1/(count + 1.5)
  xticks = range(count)
  plt = pl.figure().add_subplot(111)
  plt.set_xticks(xticks)
  plt.set_xticklabels(xticks)
  plt.set_xlabel("person giving")
  plt.set_ylabel("probability")
  plt.set_xlim(-0.5, count + 1.5)
  for i, a_probabilities in enumerate(probabilites):
    plt.bar([x + (i - (count - 1)/2) * width for x in xticks], [p/total for p in a_probabilities], width, label = i)
  plt.legend(title = "giving to", loc = "upper right")

def probability_paired(count):
  totalProbability = 0
  pairedProbability = 0
  totalCount = 0
  pairedCount = 0
  for p, entry in generate_options(count):
    if entry is not None:
      totalProbability = totalProbability + p
      totalCount = totalCount + 1
      valid = True
      for i in range(count):
        if entry[entry[i]] != i:
          valid = False
          break
      if valid:
        pairedProbability = pairedProbability + p
        pairedCount = pairedCount + 1
  return pairedProbability, totalProbability, pairedCount, totalCount

Answer 4

As @Juho-Kokkala and others have stated, the valid derangements are not equally likely and it is not stated in the problem what happens if the last person (H in the example) is the same as the last name in the bowl. I think it is understood from the question what happens when a person prior to the last one draws their own name: that person- and that person only- redraws and then replaces their name (that's equivalent to replacing their name and redrawing, but more efficient).

I suspect that in practice, the group would have the last person switch with the name chosen by the previous person (Person G in the example). That always gives a valid assignment because G did not pick H (so H, after the switch will have a name other than their own) and G will now have H after the switch. A second option to handle the case when the last person has only their own name to draw is to have all people put the names back in the bowl and redraw from the beginning. These two options give different probabilities and the valid derangements are not equally likely using either option.

The case of 8 names is not easy to analyze, but the case of 4 names is not too hard. This is the way of choosing the names: i) person A chooses until they get a name different than A. Each name {B, C, D} is equally likely. ii) person B chooses until they get a name different than B. If A chose B, then there are 3 equally likely choices; if A did not choose B, then there are 2 equally likely choices. iii) C chooses until they get a name different from C. iv) If D is the only name left, then D swaps names with C; otherwise D gets the last name remaining. That's option 1. Option 2 is identical until reaching step iv) which becomes: If D is the only name left, then replace all the names into the bowl and redraw from the beginning.

Option 1
There are these possible derangements with the probability of choosing that derangement using option 1 given next to each:
BADC $\frac{1}3\frac{1}3$
BCDA $\frac{1}3\frac{1}3\frac{1}2$
BDAC $\frac{1}3\frac{1}3$
CADB $\frac{1}3\frac{1}2\frac{1}2$
CDAB $\frac{1}3\frac{1}2\frac{1}2$
CDBA $\frac{1}3\frac{1}2\frac{1}2$
DABC $\frac{1}3\frac{1}2$
DCAB $\frac{1}3\frac{1}2\frac{1}2$
DCBA $\frac{1}3\frac{1}2\frac{1}2$

Then, the perfect pairings are BADC, CDBA, and DCBA and the probability of having a perfect pairing is $\frac{1}3\frac{1}3+\frac{1}3\frac{1}2\frac{1}2+\frac{1}3\frac{1}2\frac{1}2=\frac{5}{18}$

Option 2
First consider what happens in the first round of drawing (that is, when all the people have drawn a valid name until the last person's turn to draw). The probability of reaching the last person with only their name in the bowl is $\frac{5}{36}$. There are these possible derangements with the probability of choosing that derangement using option 2 in the first round given next to each:
BADC $\frac{1}3\frac{1}3$
BCDA $\frac{1}3\frac{1}3\frac{1}2$
BDAC $\frac{1}3\frac{1}3$
CADB $\frac{1}3\frac{1}2\frac{1}2$
CDAB $\frac{1}3\frac{1}2\frac{1}2$
CDBA $\frac{1}3\frac{1}2\frac{1}2$
DABC $\frac{1}3\frac{1}2$
DCAB $\frac{1}3\frac{1}2\frac{1}2$
DCBA $\frac{1}3\frac{1}2\frac{1}2$

Now, those probabilities add up to $\frac{31}{36}$. Therefore, accounting for the fact that the whole group will redraw in the cases where D has only their own name left to draw, those probabilities have to be divided by the sum in order to find the probability of achieving each derangement eventually (either in the first round or after re-starting as often as needed to get a valid derangement).

Finally, for option 2, the probability of having a perfect pairing is $\frac{\frac{1}9+\frac{1}{12}+\frac{1}{12}}{\frac{31}{36}}=\frac{10}{31}$

Back to the original question- the following R program enumerates all of the possible outcomes with 8 people and the probabilities using option 1 and then using option 2. The probability of a perfect pair is 0.006254409 using option 1 and 0.006876715 using option 2. This program only works for the case 8 people and would not be easy to adapt to the general case.

x=data.frame(A=as.character(rep("A",14833)),B=as.character(rep("A",14833)),C=as.character(rep("A",14833)),D=as.character(rep("A",14833)),
             E=as.character(rep("A",14833)),F1=as.character(rep("A",14833)),G=as.character(rep("A",14833)),H=as.character(rep("A",14833)),
             probinv=rep(0,14833),perfectpair=rep(F,14833),stringsAsFactors = F)
i=0
fulllist=c("A","B","C","D","E","F","G","H")
proba=7
for (a in c("B","C","D","E","F","G","H")) {
  validlistb=setdiff(fulllist,c(a,"B"))
  probb=proba*length(validlistb)
  for (b in validlistb) {
    validlistc=setdiff(fulllist,c(a,b,"C"))
    probc=probb*length(validlistc)
    for (c1 in validlistc) {
      validlistd=setdiff(fulllist,c(a,b,c1,"D"))
      probd=probc*length(validlistd)
      for (d in validlistd) {
        validliste=setdiff(fulllist,c(a,b,c1,d,"E"))
        probe=probd*length(validliste)
        for (e in validliste) {
          validlistf=setdiff(fulllist,c(a,b,c1,d,e,"F"))
          probf=probe*length(validlistf)
          for (f in validlistf) {
            validlistg=setdiff(fulllist,c(a,b,c1,d,e,f,"G"))
            if (is.element("H",validlistg)) {
              i=i+1
              x[i,1:7]=c(a,b,c1,d,e,f,"H")
              x[i,8]=setdiff(fulllist,x[i,1:7])
              x$probinv[i]=probf
          if (sum(x[i,rank(x[i,1:8])]==fulllist)==8) x$perfectpair[i]=T
            } else if (length(validlistg)==1) {
              i=i+1
              x[i,1:7]=c(a,b,c1,d,e,f,validlistg)
              x[i,8]=setdiff(fulllist,x[i,1:7])
              x$probinv[i]=probf
          if (sum(x[i,rank(x[i,1:8])]==fulllist)==8) x$perfectpair[i]=T
            } else {
              i=i+1
              x[i,1:7]=c(a,b,c1,d,e,f,validlistg[1])
              x[i,8]=setdiff(fulllist,x[i,1:7])
              x$probinv[i]=2*probf
          if (sum(x[i,rank(x[i,1:8])]==fulllist)==8) x$perfectpair[i]=T
              i=i+1
              x[i,1:7]=c(a,b,c1,d,e,f,validlistg[2])
              x[i,8]=setdiff(fulllist,x[i,1:7])
              x$probinv[i]=2*probf
          if (sum(x[i,rank(x[i,1:8])]==fulllist)==8) x$perfectpair[i]=T
            }
          }
        }
      }
    }
  }
}
sum(1/x$probinv)
sum(((1/x$probinv[x$perfectpair])))



#option 2
x=data.frame(A=as.character(rep("A",14833)),B=as.character(rep("A",14833)),C=as.character(rep("A",14833)),D=as.character(rep("A",14833)),
             E=as.character(rep("A",14833)),F1=as.character(rep("A",14833)),G=as.character(rep("A",14833)),H=as.character(rep("A",14833)),
             probinv=rep(0,14833),perfectpair=rep(F,14833),stringsAsFactors = F)
i=0
fulllist=c("A","B","C","D","E","F","G","H")
proba=7
for (a in c("B","C","D","E","F","G","H")) {
  validlistb=setdiff(fulllist,c(a,"B"))
  probb=proba*length(validlistb)
  for (b in validlistb) {
    validlistc=setdiff(fulllist,c(a,b,"C"))
    probc=probb*length(validlistc)
    for (c1 in validlistc) {
      validlistd=setdiff(fulllist,c(a,b,c1,"D"))
      probd=probc*length(validlistd)
      for (d in validlistd) {
        validliste=setdiff(fulllist,c(a,b,c1,d,"E"))
        probe=probd*length(validliste)
        for (e in validliste) {
          validlistf=setdiff(fulllist,c(a,b,c1,d,e,"F"))
          probf=probe*length(validlistf)
          for (f in validlistf) {
            validlistg=setdiff(fulllist,c(a,b,c1,d,e,f,"G"))
            probg=probf*length(validlistg)
            for (g in validlistg) {
              h=setdiff(fulllist,c(a,b,c1,d,e,f,g))
              if (h!="H") {
                i=i+1
                x[i,1:8]=c(a,b,c1,d,e,f,g,h)
                x$probinv[i]=probg
            if (sum(x[i,rank(x[i,1:8])]==fulllist)==8) x$perfectpair[i]=T
              }
            }
          }
        }
      }
    }
  }
}
sum(1/x$probinv)
sum(((1/x$probinv[x$perfectpair])/sum(1/x$probinv)))