Hats and capricious secretary: a MonteCarlo trick to compute “e”

A popular question on basic probability  appearing in quant job interviews is the classical matching problem described in amusing forms: hats in a check-room are mixed and randomly distributed to the guests; a secretary has mixed a set of letters and their envelopes. What is the probability that no hat will match the right guest? What is the probability that exactly 3 letters will reach the right destination?

The more classical formulation of the same problem is about seeking for a match among two identical decks of N cards randomly mixed. Being A_i the probability that the i^{th} card will match (regardless to other possible matchings!), we have

P(A_i) = \frac{1}{N}\,\,, P(A_iA_j) = \frac{1}{N(N-1)}\,\,, P(A_{i_1}...A_{i_{i_r}}) = \frac{(N-r)!}{N!}
The probability of at least one match is given by
P(A_1\cup A_2 \cup ... \cup A_N)

= \sum_{i=1}^{N} P(A_i)- \sum_{i<j}P(A_iA_j) +...+(-1)^{N-1}P(A_1...A_N)

= 1- \frac{1}{2!}+\frac{1}{3!}-...+(-1)^{N-1}\frac{1}{N!}
Therefore, in the limit of large N the probability approaches 1-1/e , which is a quite surprisingly result. This result can be used to compute the probability of zero matches which is simply 1 minus the probability of at least one match.

From this simple result the amusing idea to try to compute e with a MonteCarlo, just as it is usually done for \pi. Here my Mathematica code which does the job:

MonteCarloNepero[NN_, L_] := Module[{mean, check, indicatorCount},
 indicatorCount = 0;
 check = (# - RandomSample[#]) &@Table[l, {l, 1, L}];
 If[FreeQ[check, 0],
 indicatorCount = indicatorCount + 1;
 , {NN}];
 mean = indicatorCount/NN;

Where NN is the number of MonteCarlo runs and L is the length of the card deck. The convergence in the number of cards is very fast so with L=15 the result is already fair. Far more slow is the convergence in NN… For example we may get

In[2]:= MonteCarloNepero[1000, 15]
Out[2]:= 2.75482

In[3]:= E // N
Out[3]= 2.71828

Not that accurate, but fun.