Libres pensées d'un mathématicien ordinaire

This post is devoted to the classical algorithm for the simulation of a discrete law \( {P=\sum_{a\in A}p_a\delta_a} \) on an at most countable set \( {A} \). One may number the atoms by selecting a bijection \( {\sigma:A\rightarrow\{1,2,\ldots\}} \). Now if we partition the real interval \( {[0,1]} \) into blocs of respective Lebesgue measure \( {p_{\sigma^{-1}(1)},p_{\sigma^{-1}(2)}} \), etc, for instance one may use the intervals

\[ I_1=[0,p_{\sigma^{-1}(1)}[, I_2=[p_{\sigma^{-1}(1)},p_{\sigma^{-1}(1)}+p_{\sigma^{-1}(2)}[,\ldots \]

and if \( {U} \) is a uniform random variable on \( {[0,1]} \), then \( {\mathbb{P}(U\in I_{\sigma(a)})=p_a} \) for every \( {a\in A} \). The classical algorithm goes then as follows: we first generate \( {U} \) using rand, next, if \( {U\leq p_{\sigma^{-1}(1)}} \) then we output \( {\sigma^{-1}(1)} \), elseif \( {U\leq p_{\sigma^{-1}(1)}+p_{\sigma^{-1}(2)}} \), we output \( {\sigma^{-1}(2)} \), etc. The cost of this algorithm is the number of comparisons \( {N} \) needed. This number is random, and follows actually the image law of \( {P} \) by the map \( {\sigma} \). In particular, we have \( {\mathbb{E}(N)=\sum_{a\in A} \sigma(a)p_{a}} \).

Any \( {\sigma} \) minimizing the average cost \( {\mathbb{E}(N)} \) satisfies to \( {p_{\sigma^{-1}(1)}\geq p_{\sigma^{-1}(2)}\geq\cdots} \). The best average cost is obtained by labeling the atoms in decreasing probability order, and its value is the mean of this relabeled discrete law.

For the geometric distribution of arbitrary mean, and for the Poisson distribution of unit mean, the atoms are already labeled decreasingly. For finite discrete distributions, the cost of the reordering is \( {|A|\log|A|} \) with high probability (quick sort, can be done once for all).

One may ask if one can simulate the Poisson distribution of arbitrary mean \( {\lambda} \) from a generator of the Poisson distribution of unit mean. The answer is positive, and relies simply on thinning. Namely, if \( {X_1,\ldots,X_{\lceil\lambda\rceil}} \) are i.i.d. Poisson random variables with unit mean, then their sum \( {S=X_1+\cdots+X_{\lceil\lambda\rceil}} \) is Poisson with mean \( {\lceil\lambda\rceil} \), and conditional on \( {\{S=n\}} \), if \( {B_1,\ldots,B_n} \) are i.i.d. Bernoulli random variables with mean \( {\lambda/\lceil\lambda\rceil} \) then \( {B_1+\cdots+B_S} \) follows the Poisson distribution of mean \( {\lambda} \) (Poisson-Binomial mixtures are Poisson).

I like very much such kind of concrete elementary mathematics, involving basic algorithmic and probabilistic concepts. Try this with your students!

This post is devoted to an algorithm due to Aart Johannes Stam for the simulation of the uniform law on the set \( {\Pi_n} \) of partitions of \( {\{1,\ldots,n\}} \). This law gives the same probability \( {1/B_n} \) to every element of \( {\Pi_n} \), where \( {B_n=\mathrm{Card}(\Pi_n)} \). In combinatorics, \( {{(B_n)}_{n\geq1}} \) are the Bell numbers. We have \( {B_1=1} \), \( {B_2=2} \), and more generally, using the convention \( {B_0=1} \),

\[ B_{n+1}=\sum_{k=0}^n\binom{n}{k}B_k \]

(here \( {k} \) stands for the number of elements outside the bloc containing \( {n+1} \)). Now, the formal series \( {G(X)=\sum_{n=0}^\infty\frac{B_n}{n!}X^n} \) satisfies to \( {G’(X)=\exp(X)G(X)} \), and thus

\[ G(X)=\exp(\exp(X)-1). \]

We recognize the Laplace transform of the Poisson law of unit mean. In particular, the Bell numbers are the moments of this law, which gives the Dobinski formula:

\[ B_n=\frac{1}{e}\sum_{k=0}^\infty\frac{k^n}{k!}. \]

To simulate the uniform law on \( {\Pi_n} \), the idea is to associate a random color to each element of \( {\{1,\ldots,n\}} \) and to decide that elements of identical color are in the same bloc. More precisely, we first generate a random integer \( {K} \) taking the value \( {k} \) with probability \( {k^n/(k!eB_n)} \) (this is a well defined law thanks to the Dobinski formula above). Next, conditional on \( {K} \), we generate i.i.d. random variables \( {C_1,\ldots,C_n} \) according to the uniform law on \( {\{1,\ldots,K\}} \). Finally, we define the random element \( {P} \) of \( {\Pi_n} \) obtained by deciding that \( {i} \) and \( {j} \) are in the same bloc iff \( {C_i=C_j} \). It turns out that \( {P} \) follows the uniform law on \( {\Pi_n} \) since for every \( {p\in\Pi_n} \) with say \( {b} \) blocs:

\[ \begin{array}{rcl} \mathbb{P}(P=p) &=&\sum_{k=b}^\infty\mathbb{P}(P=p|K=k)\mathbb{P}(K=k)\\ &=&\sum_{k=b}^\infty\frac{k(k-1)\cdots(k-b+1)}{k^n}\frac{k^n}{k!eB_n}\\ &=&\frac{1}{B_n}. \end{array} \]

Exercice. Evaluate the complexity of this algorithm!

References.

• Stam, Generation of a random partition of a finite set by an urn model. J. Combin. Theory Ser. A 35 (1983), no. 2, 231–240;
• Knuth, The art of computer programming. Vol. 4, Fasc. 3. Generating all combinations and partitions. Addison-Wesley, 2005.

Note. I’ve learned this algorithm while writing an expository text with Yan Doumerc and Florent Malrieu on the chinese restaurant process. It is worth noting that the uniform law on \( {\Pi_n} \) does not belong to the Ewens one parameter family of laws on partitions. This is in contrast with the uniform law on permutations, which coincides with the Ewens law of unit parameter.

Note. One should not confuse finite set partitions with integer partitions, which are rather related to Young or Ferrers diagrams.

En 1965 déjà, Pierre Dac présentait sa candidature aux présidentielles, au nom du Mouvement Ondulatoire Unifié (MOU) dont la devise était : « Les temps sont durs, vive le MOU ! »

In quantum mechanics, the spin-statistics theorem states that sub-atomic particles are either bosons or fermions. Photons are bosons, while electrons are fermions. Bosons have integer spin, and the wave function of a system of two bosons is symmetric. Fermions have half integer spin, and the wave function of a system of two fermions is antisymmetric. Bosons can occupy the same quantum state, and obey to the Bose-Einstein statistics. Fermions cannot occupy the same quantum state (Pauli exclusion principle), and obey to the Fermi-Dirac statistics. These statistics can be retained using some elementary combinatorics, in the spirit of the famous Maxwell-Boltzmann statistics.

Maxwell-Boltzmann statistics. We have \( {m} \) energy levels \( {E_1,\ldots,E_m} \) and \( {n} \) distinguishable particles. If \( {n_i} \) is the number of \( {E_i} \) particles then \( {n=n_1+\cdots+n_m} \) and the average energy of configuration \( {(n_1,\ldots,n_m)} \) is \( {E=n_1E_1+\cdots+n_mE_m} \). There are

\[ \binom{n}{n_1,\ldots,n_m}=n!\prod_{i=1}^m\frac{1}{n_i!} \]

ways to obtain the configuration \( {(n_1,\ldots,n_m)} \) (multinomial coefficient: throw \( {n} \) distinguishable balls into \( {m} \) boxes). In order to get an additive measurement of the degree of freedom, we may take the logarithm. The additive degree of freedom per particle is then \( {\frac{1}{n}\log\binom{n}{n_1,\ldots,n_m}} \). Now, if \( {n\rightarrow\infty} \) with \( {n_i/n\rightarrow p_i} \) for every \( {1\leq i\leq n} \) then \( {p_1+\cdots+p_m=1} \), and thanks to the Stirling formula, the asymptotic degree of freedom per particle is given by

\[ -\sum_{i=1}^mp_i\log(p_i). \]

This quantity is known as the Boltzmann entropy. Maximizing this entropy subject to the fixed average energy constraint \( {E=\sum_{i=1}^mp_iE_i} \) gives

\[ p_i=Z^{-1}e^{-\beta E_i} \]

where the Lagrange multiplier \( {\beta>0} \) can be interpreted to be proportional to an inverse temperature, and were \( {Z=\sum_{i=1}^me^{-\beta E_i}} \) is a normalizing factor.

As pointed out by Gibbs, the reasoning above is faulty. Namely, if we consider a system formed by \( {n+n’} \) particles, then the logarithmic degree of freedom is actually not additive, violating the second principle of thermodynamics. A way to solve this paradox is to assume that the particles are undistinguishable. This makes the combinatorics trivial since for every configuration \( {(n_1,\ldots,n_m)} \) there is a unique way to throw our undistinguishable particles into the \( {m} \) boxes with \( {n_i} \) particles in each box for every \( {1\leq i\leq m} \). To escape from this triviality, we assume that each energy level \( {E_i} \) has \( {s_i} \) states, in other words that box \( {i} \) has \( {s_i} \) sub-boxes. The number of ways to obtain configuration \( {(n_1,\ldots,n_m)} \) is then

\[ \prod_{i=1}^m\frac{(n_i+s_i-1)!}{n_i!(s_i-1)!}. \]

If \( {s_1=\cdots=s_m=1} \) then this number is \( {1} \). If \( {n_i\ll s_i} \) then thanks to the Stirling formula, the number of ways to obtain configuration \( {(n_1,\ldots,n_m)} \) is

\[ \approx\prod_{i=1}^m\frac{s_i^{n_i}}{n_i!}. \]

The variational reasoning used for the distinguishable case leads then to another distribution for the energy levels, known as the Maxwell-Boltzmann statistics:

\[ p_i\approx Z^{-1}s_ie^{-\beta E_i}. \]

When \( {s_1=\cdots=s_m} \) we recover the same distribution.

Bose-Einstein statistics. We follow the reasoning used for the Maxwell-Boltzmann statistics, without the assumption \( {n_i\ll s_i} \). This leads to

\[ p_i\approx Z^{-1}\frac{s_i}{e^{\beta E_i-\alpha}-1}. \]

Fermi-Dirac statistics. Suppose this time that at most one particle can occupy a state (i.e. a sub-box). The number of ways to obtain configuration \( {(n_1,\ldots,n_m)} \) is this time

\[ \prod_{i=1}^m\frac{s_i!}{n_i!(s_i-n_i)!}, \]

leading to

\[ p_i\approx Z^{-1}\frac{s_i}{e^{\beta E_i-\alpha}+1}. \]