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.
 A: 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.
A: 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:

A: I am not sure about what I'm saying, but perhaps this could turn out to be true : 
Count the number of possibilites that gives rise to $X \ge k$. You have to choose which bins get $1$ ball, and then you have to choose which balls get in those bins. This gives you $\begin{pmatrix} n \\ k \end{pmatrix}\begin{pmatrix} n \\ k \end{pmatrix}$ possibilites. 
EDIT because of a good comment =) : You also have to choose here which ball goes in which bin, hence an extra factor of $k!$.
Now the remaining $n-k$ balls go anywhere in the $n-k$ bins left, so that gives you an extra factor of $(n-k)^{n-k}$. The total number of possibilites is just $n^n$, hence
$$
\mathbb P(X\ge k) = \frac{ \begin{pmatrix} n \\ k \end{pmatrix}\begin{pmatrix} n \\ k \end{pmatrix} k! (n-k)^{n-k} }{n^n}.
$$
To get $\mathbb P(X=k)$, just compute $\mathbb P(X \ge k) - \mathbb P(X \ge {k+1})$. 
I don't know if our formulas are asymptotically equivalent (yours and mine), but perhaps if you're more interested in this question than I am, you could try working it out. =)
Hope that helps,
P.S. After reading the comments my argument feels wrong, but I'm still going to leave it there for readers.
