A little quirk on the good old "how many ways are there to place balls into boxes" Given $n$ (unlabeled) balls, and $m$ (labeled) containers of equal size, I want to calculate what is the most probable distribution of balls, from the point of view of have many containers will remain empty. Balls are assigned randomly , with equal probabilites for each container.
So, I believe the number of configurations in which $k$ containers out of the available $m$ contain at least one ball equals
$$ C(k) = { {n -1} \choose {k -1} }{ {m} \choose {k} } $$ 
where the first term gives the number of combinations to put $n$ balls in $k$ containers with each container containing at least one ball (Put $N$ identical balls into $m$ different buckets, each bucket has at least one ball, how many ways?), and the second reflects the freedom in choosing the $k$ containers to place balls in.
It is a matter now to maximise the function $C(k)$. 
I checked that when Stirling's approximation applies, for large arguments of the factorial, the maximum tends to $n$, so filling all the boxes. In other words, with many particles, it is unlikely any box will remain empty.
But plotting $C(k)$ for smaller $n$ and $m$, one notices the maximum occurs for $k<m$, so some containers are indeed empty.
Now, the question.
I am interested in estimating when the transition occurs, i.e. for which values of $m, n$ the maximum of $C(k)$ occurs for $k < m$. I am stuck a bit as I cannot manipulate the binomial coefficients and the factorials for small values of their argument.
My attempt followed this route.  Let us assume that with given $n$ and $m$, one computes the most likely outcome is such that $k$ containers only are filled. But then, one could divide the containers in two sets, filled and empty by the process described above, each considered as a container in itself. One can view the process as a Bernoulli process, in which balls are assigned to a larger Container, union of  the selected $k$ "filled" containers) with probability proportional to $k$ (and this is a Bernoulli "success"), and the union of the $m-k$ empty containers, with probability proportional to $m-k$. The analogue of getting $k$ empty containers in the above described process is, for the Bernoulli process, to get all successes over the $n$ attempts.  For the event "all successes" to be the expected value, the following conditions has to hold
$$ n p \geq n -1/2 $$
where $p$, the chance of success on the single  draw, equals 
$$ p = \frac {m-k}{m}$$
I end up then with the condition
$$ n (\frac {m-k}{m} ) \geq n -1/2 $$
At least, I believe this confirms that when $n \to \infty$, $k \to 0$.
Any hint on how to substantiate the condition found with a direct calculation on $C(k)$ would be greatly appreciated, thanks.
 A: Using the Generalized Inclusion-Exclusion Principle
The sum of the probabilities of missing $j$ particular bins over all selections of $j$ particular bins is
$$
\overbrace{\ \ \ \binom{m}{j}\ \ \ }^{\substack{\text{the number of}\\\text{ways to choose}\\\text{$j$ particular bins}}}\overbrace{\left(\frac{m-j}m\right)^n}^{\substack{\text{the probability}\\\text{of missing $j$}\\\text{particular bins}}}\tag{1}
$$
Then, according to the Generalized Inclusion-Exclusion Principle, the probability of missing exactly $k$ bins is
$$
\begin{align}
\sum_{j=0}^m(-1)^{j-k}\binom{j}{k}\binom{m}{j}\left(\frac{m-j}m\right)^n
&=\binom{m}{k}\sum_{j=0}^m(-1)^{j-k}\binom{m-k}{m-j}\left(\frac{m-j}m\right)^n\tag{2}
\end{align}
$$
In the question, $k$ is the number of non-empty bins, so we need to substitute $k\mapsto m-k$. Since $j$ is a dummy variable, we can also substitute $j\mapsto m-j$. This gives the probability of getting exactly $k$ non-empty bins to be
$$\newcommand{\stirtwo}[2]{\left\{#1\atop#2\right\}}
\begin{align}
\binom{m}{k}\sum_{j=0}^m(-1)^{k-j}\binom{k}{j}\left(\frac{j}m\right)^n
&=\binom{m}{k}\stirtwo{n}{k}\frac{k!}{m^n}\\
&=\bbox[5px,border:2px solid #C0A000]{\stirtwo{n}{k}\frac{(m)_k}{m^n}}\tag{3}
\end{align}
$$
where $\stirtwo{n}{k}$ is a Stirling Number of the Second Kind and $(m)_k$ is a Pochhammer Symbol or Falling Factorial.
One of the "defining" relations of the Stirling Number of the Second Kind is that
$$
\sum_{k=0}^n\stirtwo{n}{k}(m)_k=m^n\tag{4}
$$
and $(4)$ guarantees that the sum of $(3)$ over all $k$ is $1$, which is comforting.
Formula $(3)$ matches user49640's formula, though I believe this approach is different.

Looking for the Maximum
Now we are looking for when $\stirtwo{n}{k}(m)_k$ is maximum for $k=m$. Since
$$
\frac{\stirtwo{n}{k}(m)_{k}}{\stirtwo{n}{k-1}(m)_{k-1}}=\frac{\stirtwo{n}{k}}{\stirtwo{n}{k-1}}(m-k+1)\tag{5}
$$
we are looking for $n$ and $m$ so that
$$
\stirtwo{n}{m}\ge\stirtwo{n}{m-1}\tag{6}
$$
According to this Wikipedia article section, the maximum $m$ satisfying $(6)$ is asymptotically
$$
m\sim\frac{n}{\log(n)}\tag{7}
$$
If some $m$ satisfies $(6)$, any smaller $m$ will.
A: This isn't a complete answer to the OP's question, but it corrects the expression given for the required probability.
Let $P(n,m,k)$ be the probability that after $n$ balls are randomly assigned to $m$ containers, exactly $k$ of them have at least one ball in them.
Using the inclusion-exclusion formula, we find
$$P(n,m,m) = 1 - m(1-1/m)^n + \binom{m}{2}(1-2/m)^n - \dots = \sum_{i=0}^{m-1} (-1)^i \binom{m}{i}(1 - i/m)^n.$$
Now, summing over all possible $k$-element subsets of the $m$ containers, we have 
$$P(n,m,k) = \binom{m}{k}(k/m)^nP(n,k,k) = \binom{m}{k} \sum_{i=0}^{k-1} (-1)^i \binom{k}{i}\left(\frac{k-i}{m}\right)^n =  \frac{m!}{m^n(m-k)!}S(n,k),$$
where $S(n,k)$ denotes a Stirling number of the second kind.
A: Suppose we have fixed value of $m$, the number of bins, and we want to find the value $n$ which maximizes the chance of exactly $k$ nonempty bins, either for a single value of $k=K$ or for all $k\le K$.  We can create a Markov process with $K+1$ states that computes these values.
Note that the labelling of the bins does not enter the computation.  If $k$ out of $m$ bins are nonempty, then drawing a ball will result either in a new bin getting a ball with probability $(m-k)/m$ or preserving the same set of bins that have at least one ball with probability $k/m$.
Assume $1 \lt K \lt m$.  Once more than $K$ bins have a ball, this will always remain true.  So we group together all states with more than $K$ nonempty bins as a single "absorbing" state.
The state after drawing one ball is deterministic: one bin is nonempty and the rest of the $m-1$ bins are empty.  So we have a state vector of length $K+1$ at step $1$:
$$ s_1 = (1,0,0,\ldots,0) $$
and the probability transition matrix looks like:
$$ M = \begin{bmatrix} 1/m & 1 - 1/m & 0 & 0 & \dots & 0 \\
0 & 2/m & 1-2/m & 0 & \dots & 0 \\
0 & 0 & 3/m & 1-3/m & \dots & 0 \\
\vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\
0 & 0 & \dots & (K-1)/m & 1-(K-1)/m & 0 \\
0 & 0 & \dots & 0 & K/m & 1-K/m \\
0 & 0 & \dots & 0 & 0 & 1 \end{bmatrix} $$
Thus the probability distribution of states at step $n$ is simply $s_1 M^{n-1}$.  The next two steps give these results:
$$ s_2 = (m^{-1},1-m^{-1},0,\ldots,0) $$
$$ s_3 = (m^{-2},3m^{-1}-3m^{-2}, 1-3m^{-1}+2m^{-2},0, \ldots,0) $$
Numerical computation to find the step $n$ that maximizes a particular entry (the $k$th entry of $s_n$ is the probability that $k$ bins are nonempty) can be scaled to integer arithmetic by extracting a factor $1/m$ from matrix $M$.  If the column vectors $e_j$, $j=1,..,K+1$ are the "standard basis vectors", then since $s_1 = e_1^T$:
$$ \Pr(k \text{ bins exactly are nonempty after } n \text{ balls}) = m^{1-n} e_1^T A^{n-1} e_k $$
where $A = mM$ is an upper triangular matrix with nonnegative integer entries.
Furthermore $A = P^{-1}D P$ is similar to its matrix $D$ of diagonal entries, so once $u^T = e_1^T P^{-1}$ and $v_k = P e_k$ are computed, the above probability is easy to find:
$$ \Pr(k \text{ bins exactly are nonempty after } n \text{ balls}) = \frac{u^T D^{n-1} v_k}{m^{n-1}} $$
