To make the problem simpler, let's consider the case where the allowed values of the share of each person is discrete, e.g., integers. Equivalently, one can also imagine partitioning the "income" axis into equally-spaced intervals and approximate the values falling into a given interval by the value at the midpoint.

Denoting the total income as $X$, the $s$-th allowed value as $x_{s}$, the total number of people as $N$, and finally, the number of people with shares of $x_{s}$ as $n_{s}$, the following conditions should be satisfied:
\begin{equation}
C_{1} (\{n_{s}\})\equiv\sum_{s} n_{s} - N = 0,
\end{equation}
and
\begin{equation}
C_{2} (\{n_{s}\})\equiv \sum_{s} n_{s} x_{s} - X = 0.
\end{equation}

Notice that many different ways to divide the share can represent the same distribution. For example, if we consider dividing \$4 between two people, giving \$3 to Alice and \$1 to Bob and vice versa both give identical distributions. As the division is random, the distribution with the maximum number of corresponding ways to divide the share has the best chance to occur.

Therefore, the goal is to maximize 
\begin{equation}
W(\{n_{s}\}) \equiv \frac{N!}{\prod_{s} n_{s}!},
\end{equation}
under the two constraints given above. The method of Lagrange multipliers is a canonical approach for this, and furthermore, one can choose to maximize $\ln W$ instead of $W$ itself, as "$\ln$" is a monotone increasing function. That is,
\begin{equation}
\frac{\partial \ln W}{\partial n_{s}} = \lambda_{1} \frac{\partial C_{1}}{\partial n_{s}} + 
\lambda_{2} \frac{\partial C_{1}}{\partial n_{s}} = \lambda_{1} + \lambda_{2} x_{s},
\end{equation}
where $\lambda_{1,2}$ are Lagrange multipliers. Notice that according to [Stirling's formula][1],
\begin{equation}
\ln n! \approx n\ln n - n,
\end{equation}
leading to
\begin{equation}
\frac{d\ln n!}{dn} \approx \ln n.
\end{equation}
Thus,
\begin{equation}
\frac{\partial \ln W}{\partial n_{s}} \approx -\ln n_{s}.
\end{equation}
It then follows that
\begin{equation}
n_{s} \approx \exp\big(-\lambda_{1} - \lambda_{2} x_{s}\big),
\end{equation}
which is an exponential distribution.

That this is really a maximum, rather than a minimum or a saddle point, can be seen from the Hessian of $\ln W - \lambda_{1} C_{1} - \lambda_{2} C_{2}$. Because $C_{1,2}$ are linear in $n_{s}$, it is the same as that of $\ln W$:
\begin{equation}
\frac{\partial^{2} \ln W}{\partial n_{s}^{2}} = -\frac{1}{n_{s}} < 0,
\end{equation}
and
\begin{equation}
\frac{\partial^{2} \ln W}{\partial n_{s}\partial n_{r}} = 0 \quad (s\neq r).
\end{equation}
Hence the Hessian is concave, and what we have found is indeed a maximum.

Notice that $W(\{n_{s}\})$ is in a sense the distribution of distributions. For the distributions we typically observe to be something close to the most probable one, $W(\{n_{s}\})$ should be narrow enough. It is seen from the Hessian that this condition amounts to $n_{s}\gg 1$. Therefore, to actually observe the most probable (exponential) distribution, the partition in the money interval should not be too fine. Even with this, the exponential distribution is destined to become inaccurate towards the tail.

**Note:** This is exactly how physicists understand the [Boltzmann distribution][2] in statistical mechanics.


  [1]: http://en.wikipedia.org/wiki/Stirling's_approximation
  [2]: http://en.wikipedia.org/wiki/Boltzmann_distribution