How many surjections are there from a set of size n? It is well-known that the number of surjections from a set of size n to a set of size m is quite a bit harder to calculate than the number of functions or the number of injections. (Of course, for surjections I assume that n is at least m and for injections that it is at most m.) It is also well-known that one can get a formula for the number of surjections using inclusion-exclusion, applied to the sets $X_1,...,X_m$, where for each $i$ the set $X_i$ is defined to be the set of functions that never take the value $i$. This gives rise to the following expression: $m^n-\binom m1(m-1)^n+\binom m2(m-2)^n-\binom m3(m-3)^n+\dots$.
Let us call this number $S(n,m)$. I'm wondering if anyone can tell me about the asymptotics of $S(n,m)$. A particular question I have is this: for (approximately) what value of $m$ is $S(n,m)$ maximized? It is a little exercise to check that there are more surjections to a set of size $n-1$ than there are to a set of size $n$. (To do it, one calculates $S(n,n-1)$ by exploiting the fact that every surjection must hit exactly one number twice and all the others once.) So the maximum is not attained at $m=1$ or $m=n$. 
I'm assuming this is known, but a search on the web just seems to lead me to the exact formula. A reference would be great. A proof, or proof sketch, would be even better.
Update. I should have said that my real reason for being interested in the value of m for which S(n,m) is maximized (to use the notation of this post) or m!S(n,m) is maximized (to use the more conventional notation where S(n,m) stands for a Stirling number of the second kind) is that what I care about is the rough size of the sum. The sum is big enough that I think I'm probably not too concerned about a factor of n, so I was prepared to estimate the sum as lying between the maximum and n times the maximum.
 A: If f is an arbitrary surjection from N onto M, then we can think of f as partitioning N into m different groups, each group of inputs representing the same output point in M. The Stirling Numbers of the second kind count how many ways to partition an N element set into m groups. But this undercounts it, because any permutation of those m groups defines a different surjection but gets counted the same. There are m! such permutations, so our total number of surjections is
m! S(n,m)
To look at the maximum values, define a sequence S_n = n - M_n where M_n is the m that attains maximum value for a given n - in other words, S_n is the "distance from the right edge" for the maximum value. Computer-generated tables suggest that this function is constant for 3-4 values of n before increasing by 1. If this is true, then the m coordinate that maximizes m! S(n,m) is bounded by n - ceil(n/3) - 1 and n - floor(n/4) + 1.
I have no proof of the above, but it gives you a conjecture to work with in the meantime.
A: I don't have a precise reference for your problem (given $n$ find "the most surjected" $m$); waiting for a precise one, I can say that I think the standard starting point  should be as follows. To avoid confusion I modify slightly your notation  for the surjections from an $n$ elements set to an $m$ elements set into $\mathrm{Sur}(n,m).$ One has the generating function (coming e.g. from the analogous g.f. for Stirling numbers of second kind)
$$(e^x-1)^m\,=\sum_{n\ge m}\ \mathrm{Sur}(n,m)\ \frac{x^n}{n!},$$
whence by the Cauchy formula with a simple integration contour around 0 ,
$$\frac{\mathrm{Sur}(n,m)}{n!}={1 \over 2\pi i} \oint \frac{(e^z-1)^m}{z^{n+1}}dz$$
For a circular path $re^{it}$ we find
$$\frac{\mathrm{Sur}(n,m)}{n!}={1 \over 2\pi } \int_{-\pi}^{\pi}\left(\exp(re^{it})-1\right)^m e^{-int}  dt\\ .$$
This holds for any number $r>0$, and the most convenient one should be chosen according to the stationary phase method; here a change of variable followed by dominated convergence may possibly give a convergent integral, producing an asymptotics: this is e.g. how one can derive the Stirling asymptotics for n!.
In your case, the problem is: for a given $n$ (large) maximize the integral in $m$, and give asymptotic expansions for the maximal $m$ (the first order should be $\lambda n + O(1)$ with $ 2/3\leq \lambda\leq  3/4 $ according to Michael Burge's exploration). This seems to be tractable; for the moment I leave this few hints hoping they are useful, but I'm very curious to see the final answer.
A: It seems to be the case that the polynomial $P_n(x) =\sum_{m=1}^n
m!S(n,m)x^m$ has only real zeros. (I know it is true that $\sum_{m=1}^n
S(n,m)x^m$ has only real zeros.) If this is true, then the value of $m$ 
maximizing $m!S(n,m)$  is within 1 of $P'_n(1)/P_n(1)$ by a theorem of 
J. N.  Darroch, Ann. Math. Stat. 35 (1964), 1317-1321. See also 
J. Pitman, J. Combinatorial Theory, Ser. A 77 (1997), 279-303. 
By standard combinatorics
  $$ \sum_{n\geq 0} P_n(x) \frac{t^n}{n!} = \frac{1}{1-x(e^t-1)}. $$
Hence 
  $$ \sum_{n\geq 0} P_n(1)\frac{t^n}{n!} = \frac{1}{2-e^t} $$
  $$ \sum_{n\geq 0} P'_n(1)\frac{t^n}{n!} = \frac{e^t-1}{(2-e^t)^2}. $$
Since these functions are meromorphic with smallest singularity at $t=\log 2$, 
it is routine to work out the asymptotics, though I have not bothered to 
do this.
Update. It is indeed true that $P_n(x)$ has real zeros. This is because
$(x-1)^nP_n(1/(x-1))=A_n(x)/x$, where $A_n(x)$ is an Eulerian polynomial. It 
is known that $A_n(x)$ has only real zeros, and the operation $P_n(x) 
\to (x-1)^nP_n(1/(x-1))$ leaves invariant the property of having real
zeros.
A: Many people may be interested in the asymptotics for $n=cm$ where $c$ is constant (say $c=2$).
That is, how likely is a function from $2m$ to $m$ to be onto?
I couldn't dig the answer out from some of the sources and answers here, but here is a way that seems okay.
Bender (Central and local limit theorems applied to asymptotics enumeration) shows
$$\begin{eqnarray*}{n\brace m}&\sim&\frac{n!e^{-\alpha m}}{m!\rho^{n-1}(1+e^\alpha)\sigma\sqrt{2\pi n}},\\
\frac{n}m &=& (1+e^\alpha)\ln(1+e^{-\alpha}),\\
\rho&=&\ln(1+e^{-\alpha}),\\
\sigma^2&=&\left(\frac{m}n\right)^2[1-e^\alpha\ln(1+e^{-\alpha})].\end{eqnarray*}$$
Notice that for constant $n/m$, all of $\alpha$, $\rho$, $\sigma$ are constants.
For $c=2$, we find $\alpha=-1.366$
and then $\rho=1.59$
and $$\sigma^2=\frac14 [1-e^\alpha\ln(1+e^{-\alpha})]=\frac14 [1-e^{-1.366}(1.59)]=.149$$
Thus the probability that our function from $cm$ to $m$ is onto is
$$\Pr(\text{onto})=\frac1{m^n}m!{n\brace m}=\frac1{m^n}m!\frac{n!e^{-\alpha m}}{m!\rho^{n-1}(1+e^\alpha)\sigma\sqrt{2\pi n}}$$
$$=\frac1{m^n}\frac{n!e^{-\alpha m}}{\rho^{n-1}(1+e^\alpha)\sigma\sqrt{2\pi n}}$$
Does it go to 0? Well, $\rho=1.59$ and $e^{-\alpha}=3.92$, so up to polynomial factors we have
$$\frac1{m^n}\frac{n!e^{-\alpha m}}{\rho^{n}}\approx\frac{n!(3.92^m)}{(1.59)^n(n/2)^n}$$
$$\approx \frac{(3.92^m)}{e^{n}(1.59)^n(1/2)^n}=\left(\frac{2^2\times 3.92}{1.59^2\times e^2}\right)^m=0.839^m$$
So phew... it goes to 0, but not as fast as for the case $n=m$ which gives $(1/e)^m$.
A: Richard's answer is short, slick, and complete, but I wanted to mention here that there is also a "real variable" approach that is consistent with that answer; it gives weaker bounds at the end, but also tells a bit more about the structure of the "typical" surjection.  I'll write the argument in a somewhat informal "physicist" style, but I think it can be made rigorous without significant effort.
Tim's function $Sur(n,m) = m! S(n,m)$ obeys the easily verified recurrence $Sur(n,m) = m ( Sur(n-1,m) + Sur(n-1,m-1) )$, which on expansion becomes
$Sur(n,m) = \sum m_1 ... m_n = \sum \exp( \sum_{j=1}^n \log m_j )$
where the sum is over all paths $1=m_1 \leq m_2 \leq \ldots \leq m_n = m$ in which each $m_{i+1}$ is equal to either $m_i$ or $m_i+1$; one can interpret $m_i$ as being the size of the image of the first $i$ elements of $\{1,\ldots,n\}$.  If we make the ansatz $m_j \approx n f(j/n)$ for some nice function $f: [0,1] \to {\bf R}^+$ with $f(0)=0$ and $0 \leq f'(t) \leq 1$ for all $t$, and use standard entropy calculations (Stirling's formula and Riemann sums, really), we obtain a contribution to $Sur(n,m)$ of the form
$\exp( n \int_0^1 \log(n f(t))\ dt + n \int_0^1 h(f'(t))\ dt + o(n) )$  (*)
where $h$ is the entropy function $h(\theta) := -\theta \log \theta - (1-\theta) \log (1-\theta)$.  So, heuristically at least, the optimal profile comes from maximising the functional
$\int_0^1 \log(f(t)) + h(f'(t))\ dt$
subject to the boundary condition $f(0)=0$.  (The fact that $h$ is concave will make this maximisation problem nice and elliptic, which makes it very likely that these heuristic arguments can be made rigorous.)  The Euler-Lagrange equation for this problem is
$-\frac{f''}{f'(1-f')} = \frac{1}{f}$
while the free boundary at $t=1$ gives us the additional Neumann boundary condition  $f'(1)=1/2$.  The translation invariance of the Lagrangian gives rise to a conserved quantity; indeed, multiplying the Euler-Lagrange equation by $f'$ and integrating one gets
$\log(1-f') = \log f + C$
which is easily solved as
$f = \frac{1}{A} (1 - B e^{-At} )$
for some constants A, B.  The Dirichlet boundary condition $f(0)=0$ gives $B=1$; the Neumann boundary condition $f'(1)=1/2$ gives $A=\log 2$, thus
$f(t) = (1 - 2^{-t}) / \log 2$.
In particular $f(1)=1/(2 \log 2)$, which matches Richard's answer that the maximum occurs when $m/n \approx 1/(2 \log 2)$.  To match up with the asymptotic for $Sur(n,m)$ in Richard's answer (up to an error of $\exp(o(n))$, I need to have
$\int_0^1 \log f(t) + h(f'(t))\ dt = - 1 - \log \log 2.$
And happily, this turns out to be the case (after a mildly tedious computation.)
This calculation reveals more about the structure of a "typical" surjection from n elements to m elements for m free, other than that $m/n \approx 1/(2 \log 2)$; it shows that for any $0 < t < 1$, the image of the first $tn$ elements has cardinality about $f(t) n$.  If one fixes $m$ rather than lets it be free, then one has a similar description of the surjection but one needs to adjust the A parameter (it has to solve the transcendental equation $(1-e^{-A})/A = m/n$).
With a bit more effort, this type of computation should also reveal the typical distribution of the preimages of the surjection, and suggest a random process that generates something that is within o(n) edits of a random surjection.  
It's also interesting to note that the answer $m/n \approx 1/(2\log 2) = 0.72134\ldots$ fits extremely well with Kevin's numerical computation $f(1000)=722$, so we now have several independent confirmations that this is the correct answer...
A: This looks like the Stirling numbers of the second kind (up to the $m!$ factor). 
This and this papers are specifically devoted to the maximal Striling numbers.
It seems that for large $n$ the relevant asymptotic expansion is
$$k! S(n,k)= (e^r-1)^k \frac{n!}{r^n}(2\pi k B)^{-1/2}\left(1-\frac{6r^2\theta^2 +6r\theta+1}{12re^r}+O(n^{-2})\right),$$
where
$$e^r-1=k+\theta,\quad \theta=O(1),$$
$$B=\frac{re^{2r}-(r^2+r)e^r}{(e^r-1)^2}.$$ 
A: I found this paper of Temme (available here) that gives an explicit but somewhat complicated asymptotic for the Stirling number S(n,m) of the second kind, by the methods alluded to in previous answers (generating functions -> contour integral -> steepest descent)
Here's the asymptotic (as copied from that paper).  One first sets
$t_0 := \frac{n-m}{m}$
and finds the positive real number $x_0$ solving the transcendental equation
$\frac{1-e^{-x_0}}{x_0} = \frac{m}{n}$
(one has the asymptotics $x_0 \approx 2(1-m/n)$ when $m/n$ is close to 1, and $x_0 \approx n/m$ when $m/n$ is close to zero.)  One then defines
$A := \phi(x_0) - m t_0 + (n-m) t_0$
where
$\phi(x) := - n \ln x + m \ln(e^x - 1).$
(Note: $x_0$ is the stationary point of $\phi(x)$.)  One has an integral representation
$S(n,m) = \frac{n!}{m!} \frac{1}{2\pi i} \int e^{\phi(x)} \frac{dx}{x}$
where the integral is a small contour around the origin.  The saddle point method then gives
$S(n,m) = (1+o(1)) e^A m^{n-m} f(t_0) \binom{n}{m}$
where
$f(t_0) := \sqrt{\frac{t_0}{(1+t_0)(x_0-t_0)}}$
and o(1) goes to zero as $n \to \infty$ (uniformly in m, I believe).  
In principle, one can now approximate $m! S(n,m)$ to within o(1) and compute its maximum in finite time, but this seems somewhat tedious.  It does seem though that the maximum is attained when $m/n = c+o(1)$ for some explicit constant $0 < c < 1$.
EDIT: Actually, it's clear that the maximum is going to be obtained in the range $n/e \leq m \leq n$ asymptotically, because $m! S(n,m)$ equals $n! \approx (n/e)^n$ when $m=n$, and on the other hand we have the trivial upper bound $m! S(n,m) \leq m^n$.  Among other things, this makes $x_0$ and $t_0$ bounded, and so the f(t_0) term is also bounded and not of major importance to the asymptotics.  The other terms however are still exponential in n...
EDIT: There is also the identity
$\sum_{k=1}^n (k-1)! S(n,k) = (-1)^n Li_{1-n}(2)$
where $Li_s$ is the polylogarithm function.  So, up to a factor of n, the question is the same as that of obtaining an asymptotic for $Li_{1-n}(2)$ as $n \to -\infty$.  This seems quite doable (presumably from yet another contour integration and steepest descent method) but a quick search of the extant asymptotics didn't give this immediately.
