With $n$ balls and $n$ bins, what is the probability that exactly $k$ bins have exactly $1$ ball?

I've got a balls and bins problem. Suppose I throw $n$ balls uniformly at random into $n$ bins. What is the probability that exactly $k$ bins end up with exactly $1$ ball?

I know this seems a classical problem and may look "simple" or "naive," but I've worked days on it and still can't get the answer.

However, I think I do have a good approximation for it. Namely, let $X$ denote the number of such bins. Then

$$ Pr(X=k) \approx \binom{n}{k}\left(\frac{1}{e}\right)^{k}\left(1-\frac{1}{e}\right)^{n-k} $$

where $1/e$ is an approximation for $(1-1/n)^{n-1}$.

This approximation works great when $n$ is big and poorly when $n$ is small (like $n<5$).

Anyway, I'm looking for an exact expression. Anyone have an idea?

PS: I've written a simple simulation in C++; you can check your answer with it first: Simulation Code Here.

Exact formula can be obtained using de Moivre's formula for occurrence of exactly $k$ exchangeable events:

$$ p_n(k) = \binom{n}{k} \sum_{i=k}^n (-1)^{i-k} \binom{n-k}{i-k} \mathbb{P}(A_1 \ldots A_i) $$

Here $A_1$ is the event that the first bin contains $1$ ball, $A_1 A_2$ is the event that first two bins each contain 1 ball and so on.

$$ \begin{eqnarray} \mathbb{P}(A_1) &=& \binom{n}{1} \frac{1}{n} \left( 1 - \frac{1}{n} \right)^{n-1} = \left( 1 - \frac{1}{n} \right)^{n-1} \\ \mathbb{P}(A_1 A_2) &=& \binom{n}{1,1,n-2} \frac{1}{n^2} \left( 1- \frac{2}{n} \right)^{n-2} = \frac{n-1}{n} \left( 1- \frac{2}{n} \right)^{n-2} \\ \mathbb{P}(A_1 \ldots A_i) &=& \binom{n}{i} i! \frac{1}{n^i} \left( 1 - \frac{i}{n} \right)^{n-i} \end{eqnarray} $$

Hence the exact result you seek is: $$ p_n(k) = \binom{n}{k} \sum_{i=k}^n (-1)^{i-k} \binom{n-k}{i-k} \binom{n}{i} i! \frac{1}{n^i} \left( 1 - \frac{i}{n} \right)^{n-i} $$

This agrees with the simulations.

Let $S_i$ be the set of outcomes with $1$ ball in bin $i$. Let $N_j$ be the number of outcomes in the intersections of $j$ of the $S_i$; e.g. $N_3=\sum_{i<j<k}|S_i\cap S_j\cap S_k|$. There are $\binom{n}{j}$ choices of the $S_i$ to intersect, for each choice of $S_i$, there are $\binom{n}{j}j!$ choices and orders of balls to put into those $j$ bins, and $(n-j)^{n-j}$ ways to arrange the other $n-j$ balls in the other $n-j$ bins. Thus, $$ N_j=\binom{n}{j}\binom{n}{j}j!(n-j)^{n-j} $$ Since there are $n^n$ possible outcomes, to compute the probability of getting exactly $k$ bins with $1$ ball, use the Generalized Inclusion-Exclusion Principle: $$ \begin{align} \frac{1}{n^n}\sum_{j=k}^n(-1)^{j-k}\binom{j}{k}N_j &=\frac{1}{n^n}\sum_{j=k}^n(-1)^{j-k}\binom{j}{k}\binom{n}{j}\binom{n}{j}j!(n-j)^{n-j}\tag{1}\\ &=\frac{1}{n^n}\sum_{j=k}^n(-1)^{j-k}\binom{n}{k}\binom{n-k}{n-j}\frac{n!}{(n-j)!}(n-j)^{n-j}\\ &=\binom{n}{k}\frac{n!}{n^n}\sum_{j=k}^n(-1)^{j-k}\binom{n-k}{n-j}\frac{(n-j)^{n-j}}{(n-j)!}\\ &=\binom{n}{k}\frac{n!}{n^n}\sum_{j=0}^{n-k}(-1)^{n-k-j}\binom{n-k}{j}\frac{j^j}{j!} \end{align} $$ Appendix:

To verify that the probabilities for $k=0,1,\dots,n$ sum to $1$, $(1)$ can be summed fairly easily in $k$: $$ \begin{align} &\sum_{k=0}^n\frac{1}{n^n}\sum_{j=k}^n(-1)^{j-k}\binom{j}{k}\binom{n}{j}\binom{n}{j}j!(n-j)^{n-j}\\ &=\frac{1}{n^n}\sum_{j=0}^n\sum_{k=0}^j(-1)^{j-k}\binom{j}{k}\binom{n}{j}\binom{n}{j}j!(n-j)^{n-j}\\ &=\frac{1}{n^n}\sum_{j=0}^n(-1)^{j}0^j\binom{n}{j}\binom{n}{j}j!(n-j)^{n-j}\\ &=\frac{1}{n^n}(-1)^00^0\binom{n}{0}\binom{n}{0}0!(n-0)^{n-0}\\ &=1 \end{align} $$ Mathematica:

Here is the plot for $80$ bins, which matches Sasha's plot: