# Repeated draws from multinomial distribution

(This is a cross-post from Math StackExchange https://math.stackexchange.com/questions/609641/multinomial-distribution-sum-of-squared-probabilities)

Let $\vec X = (X_1, \dots, X_k)$ be a draw from a fair multinomial distribution with $n$ trials, i.e. $P(X_1 = x_1, \dots, X_k = x_k) = \binom{n}{x_1, \dots, x_k} k^{-n}$

Let $\vec Y$ be an independent draw from the same distribution, i.e. $\vec X$ and $\vec Y$ follow the same law.

My questions is this: what is the probability that $\vec X = \vec Y$?

It appears that this probability, for fixed $k$, is roughly $(c_k + o(1)) n^{-(k-1)/2}$. This makes sense; the first $k-1$ coordinates vary on the order of $\sqrt{n}$ and are "independent-ish". (The last coordinate is of course fixed by the first $k-1$). Can one compute $c_k$?

This is essentially a product of local CLTs (since one can uncover the multinomial distribution variable by variable: decide first how many of the $n$ variables equal "1" with probability of success $1/k$ for each one and thus variance $(n/k)(1-1/k)$; Then, of the roughly $n(k-1)/k$ remaining variables, decide how many are "2" with probability of success $1/(k-1)$ each, and variance $(n/k)(1-1/(k-1))$, and so on).

** Edit: As pointed out by Will Nelson, there is a factor of $1/\sqrt{2}$ in each step that comes from the fact that we are using the local CLT not at the mean but rather at the location of the first multinomial $X$.

Thus, the probability is asymptotic to $$\left(\frac{k}{4\pi n}\right)^{(k-1)/2}\cdot \frac{1}{\sqrt{\prod_{j=2}^k (1-j^{-1})}}$$

For $k=2$, this is just coming from $$\left(\frac{1}{\sqrt{2\pi n\cdot (1/4)}}\right)^2\int e^{-2\cdot (x^2/(2\cdot (1/4))}dx=1/\sqrt{\pi n}$$

• In the answer I posted, the $k=2$ value is computed to be $1/\sqrt{\pi}$, which differs from this answer by a factor of $\sqrt{2}$. Do you see how this might be? I wrote a quick simulation which seems to agree with the value $1/\sqrt{\pi}$. – Will Nelson Dec 28 '13 at 13:11
• Thanks for pointing out the discrepancy. Yes, I see where this factor comes from. My computation has a flaw in it, in that I approximated too crudely the term coming from the exponential in the Gaussian distribution (in the local CLT). When one takes that into account one recovers (for $k=2$) a factor of $1/\sqrt{2}$, which coincides with your answer. I will update the answer to reflect this. – ofer zeitouni Dec 28 '13 at 21:22

I roughly outline a computation for $k=2$ below. Assuming I have not made a mistake for $k=2$, generalization to $k>2$ will likely be straightforward. If there seems to be demand, I plan on posting a proof of at least the $k=2$ case when I get a chance.

Direct computer simulation suggests $c_2 \approx 0.564$, which is roughly $\frac{1}{\sqrt{\pi}}$.

Here is the computation (again, this isn't a proof, but it's in the right direction): \begin{eqnarray} \sqrt{n}\sum_{i=0}^n \frac{1}{2^{2n}} \left(\begin{array}{c} n \\ i\end{array}\right)^2 &=& \frac{1}{2\pi\sqrt{n}}\left( \sum_{i=0\ ; \ \ x\equiv i/n}^n \frac{1}{x (1-x)} e^{-2 n (x\log x + (1-x) \log(1-x) + \log 2)} + O(n^{-1})\right) \\ &=& \frac{1}{2\pi\sqrt{n}}\left( \sum_{i=0\ ; \ \ x\equiv i/n}^n 4 e^{-4 n (x-\frac{1}{2})^2} + O(n^{-1}) \right) \\ &\to& \frac{4}{2\pi} \int_{-\infty}^{\infty} e^{-4 x^2} \ dx \\ &=& \frac{1}{\sqrt{\pi}}. \end{eqnarray} The first line uses Stirling's approximation. The second uses a Taylor expansion at $x=\frac{1}{2}$. The third utilizes a straightforward limit of an approximation of the Riemann integral given on the third line.

For $k>2$, the same direct calculation approach yields: $$c_k = \frac{k^{\frac{k}{2}}}{(4\pi)^{\frac{k-1}{2}}}.$$ Not surprisingly, this is the same result obtained by @oferzeitouni (though simplified). The local CLT approach may ultimately be simpler, but the straight analytic proof from Stirling's formula and a Taylor expansion is not difficult.

• For $k=2$, you can use Stirling's approximation on ${2n \choose n} = \sum_{i=0}^n {n \choose i}^2$. – Douglas Zare Dec 28 '13 at 13:23
• To make this rigorous for higher dimensions, I think you want to use a local central limit theorem for lattice random variables. – Douglas Zare Dec 28 '13 at 16:47
• @DouglasZane Nice observation about the $\left(\begin{array}{c}2n\\ n\end{array}\right)$ identity. It yields a very simple proof for the $k=2$ case. – Will Nelson Dec 29 '13 at 7:50
• Using the exact values calculated in my answer, I can confirm that the general value of $c_k$ you give is correct to high precision. This is for fixed $k$ and $n\to\infty$. If $k$ and $n$ are related, such as $k=n$, these asymptotics are not correct. – Brendan McKay Dec 29 '13 at 8:47
• @BrendanMcKay That's encouraging, thanks. The original question asks to compute $c_k$ "for fixed $k$", so these asymptotics are clearly what is asked for. – Will Nelson Dec 29 '13 at 9:03

Write the probability as $k^{-2n}Q(k,n)$, where $Q(k,n)$ is the sum of the squares of the multinomial coefficients. With the help of OEIS we find that these have been studied before:

The first few links give single summations and asymptotic formulas. The last link says that $Q(k,n)$ is something called $W_k(2n)$ for which several references are given, and a generating function which I am guess-generalising here: $$\sum_{n=0}^\infty \frac{Q(k,n)}{n!^2}x^n = \left( \sum_{j=0}^\infty \frac{x^j}{j!^2} \right)^k.$$ I didn't prove this but I bet it is proved in one of the references given in the last link. (Actually I think it is elementary: equate coefficients on each side with the original formulation of the problem in mind.) It looks like the ants' pants for deriving rigorous asymptotics. The thing inside the large parens is of course a Bessel function.

This also gives a recurrence: $$Q(k,n) = n!^2 \sum_{j=0}^n \frac{Q(k-1,n-j)}{j!^2 (n-j)!^2},$$ starting with $Q(0,n)=\delta_{0,n}$.