How can I analytically prove that randomly dividing an amount results in an exponential distribution (of e.g. income and wealth)?

Question

In this current article in SCIENCE the following is being proposed:

Suppose you randomly divide 500 million in income among 10,000 people. There's only one way to give everyone an equal, 50,000 share. So if you're doling out earnings randomly, equality is extremely unlikely. But there are countless ways to give a few people a lot of cash and many people a little or nothing. In fact, given all the ways you could divvy out income, most of them produce an exponential distribution of income.

I have done this with the following R code which seems to reaffirm the result:

library(MASS)

w <- 500000000 #wealth
p <- 10000 #people

d <- diff(c(0,sort(runif(p-1,max=w)),w)) #wealth-distribution
h <- hist(d, col="red", main="Exponential decline", freq = FALSE, breaks = 45, xlim = c(0, quantile(d, 0.99)))

fit <- fitdistr(d,"exponential")
curve(dexp(x, rate = fit$estimate), col = "black", type="p", pch=16, add = TRUE)

My question
How can I analytically prove that the resulting distribution is indeed exponential?

Addendum
Thank you for your answers and comments. I have thought about the problem and came up with the following intuitive reasoning. Basically the following happens (Beware: oversimplification ahead): You kind of go along the amount and toss a (biased) coin. Every time you get e.g. heads you divide the amount. You distribute the resulting partitions. In the discrete case the coin tossing follows a binomial distribution, the partitions are geometrically distributed. The continuous analogues are the poisson distribution and the exponential distribution respectively! (By the same reasoning it also becomes intuitively clear why the geometrical and the exponential distribution have the property of memorylessness - because the coin doesn't have a memory either).

Answer 1

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 approximating all values falling into a given interval by 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 considered dividing \$4 between two people, giving \$3 to Alice and \$1 to Bob and vice versa would 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.

To obtain such a distribution, one has 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. Furthermore, one can choose to work with $\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, \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. One can obtain the values of Lagrange multipliers using the constraints. From the first constraint, \begin{equation} \begin{split} N &= \sum_{s} n_{s} \approx \sum_{s} \exp\big(-\lambda_{1} - \lambda_{2} x_{s}\big)\\ &\approx \frac{1}{\Delta x} \int_{0}^{\infty} \exp\big(-\lambda_{1} - \lambda_{2} x\big) \,\,dx\\ &=\frac{1}{\lambda_{2}\Delta x}\exp\big(-\lambda_{1}\big), \end{split} \end{equation} where $\Delta x$ is the spacing between allowed values. Similarly, \begin{equation} \begin{split} X &= \sum_{s} n_{s}x_{s} \approx \sum_{s} x_{s}\,\exp\big(-\lambda_{1} - \lambda_{2} x_{s}\big)\\ &\approx \frac{1}{\Delta x} \int_{0}^{\infty} x\,\exp\big(-\lambda_{1} - \lambda_{2} x\big) \,\,dx\\ &=\frac{1}{\lambda_{2}^{2}\,\Delta x}\exp\big(-\lambda_{1}\big). \end{split} \end{equation} Therefore, we have \begin{equation} \exp\big(-\lambda_{1}\big) = \frac{N^{2} \Delta x}{X}, \end{equation} and \begin{equation} \lambda_{2} = \frac{N}{X}. \end{equation} 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.

The function $W(\{n_{s}\})$ is really the distribution of distributions. For distributions we typically observe to be 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$. (It is also the condition that Stirling's formula is reliable.) Therefore, to actually see the exponential distribution, partitions in the income axis (corresponding to bins in OP's histogram) should be wide enough so that number of people in a partition is much greater than unity. Towards the tail, where $n_{s}$ tends to zero, this condition is always destined to fail.

Note: This is exactly how physicists understand the Boltzmann distribution in statistical mechanics. The exponential distribution is essentially exact for this case, as we consider $N\sim 10^{23}$.

Answer 2

In fact you can prove it's not actually exponential, almost trivially:

Compute the probability that a given share is greater than $500$ million. Compare with the probability that an exponential random variable is greater than $500$ million.

However, it's not too hard to see that for your uniform-gap example that it should be close to exponential.

Consider a Poisson process - where events occur at random over along some dimension. The number of events per unit of the interval has a Poisson distribution, and the gap between events is exponential.

If you take a fixed interval then the events in a Poisson process that fall within it are uniformly distributed in the interval. See here.

[However, note that because the interval is finite, you simply can't observe larger gaps than the interval length, and gaps nearly that large will be unlikely (consider, for example, in a unit interval - if you see gaps of 0.04 and 0.01, the next gap you see can't be bigger than 0.95).]

So apart from the effect of restricting attention to a fixed interval on the distribution of the gaps (which will reduce for large $n$, the number of points in the interval), you would expect those gaps to be exponentially distributed.

Now in your code, you're dividing the unit interval by placing uniforms and then finding the gaps in successive order statistics. Here the unit interval is not time or space but represents a dimension of money (imagine the money as 50000 million cents laid out end to end, and call the distance they cover the unit interval; except here we can have fractions of a cent); we lay down $n$ marks, and that divides the interval into $n+1$ "shares". Because of the connection between the Poisson process and uniform points in an interval, the gaps in the order statistics of a uniform will tend to look exponential, as long as $n$ is not too small.

More specifically, any gap that starts in the interval placed over the Poisson process has a chance to be "censored" (effectively, cut shorter than it would otherwise have been) by running into the end of the interval.

$\hspace{1cm}$

Longer gaps are more likely to do that than shorter ones, and more gaps in the interval means the average gap length must go down -- more short gaps. This tendency to be 'cut off' will tend to affect the distribution of longer gaps more than short ones (and there's no chance any gap limited to the interval will exceed the length of the interval -- so the distribution of gap size should decrease smoothly to zero at the size of the whole interval).

In the diagram, a longish interval at the end has been cut shorter, and a relatively shorter interval at the start is also shorter. These effects bias us away from exponentiality.

(The actual distribution of the gaps between $n$ uniform order statistics is Beta(1,n). )

So we should see the distribution at large $n$ look exponential in the small values, and then less exponential at the larger values, since the density at its largest values will drop off more quickly.

Here's a simulation of the distribution of gaps for n=2:

Not very exponential.

But for n=20, it starts to look pretty close; in fact as $n$ grows large it will be well approximated by an exponential with mean $\frac{1}{n+1}$.

If that was actually exponential with mean 1/21, then $\exp(-21x)$ would be uniform... but we can see it isn't, quite:

The non-uniformity in the low values there corresponds to large values of the gaps -- which we'd expect from teh above discussion, because the effect of the "cutting off" the Poisson process to a finite interval means we don't see the largest gaps. But as you take more and more values, that goes further out into the tail, and so the result starts to look more nearly uniform. At $n=10000$, the equivalent display would be harder to distinguish from uniform - the gaps (representing shares of the money) should be very close to exponentially distributed except at the very unlikely, very very largest values.

Answer 3

Let's suppose the money is infinitely divisible so we can deal with real numbers rather than integers.

Then the uniform distribution of $t=500000000$ partitioned across $n=10000$ individuals will give a marginal density for each individual $$p(x)=\frac{n-1}{t}\left(1-\frac{x}{t}\right)^{n-2}$$ for $0 \le x \le t$, and a marginal cumulative probability for each individual of $$P(X \le x)=1 - \left(1-\frac{x}{t}\right)^{n-1}.$$

If you want to apply this then use the marginal distribution to allocate a random amount $X$ to any of the individuals, then reduce $t$ to $t-X$ and $n$ to $n-1$ and repeat. Note that when $n=2$, this would give each individual a uniform marginal distribution across the remaining amount, much as one might expect; when $n=1$ you give all the remaining money to the single remaining person.

These expressions are polynomial rather than exponential, but for large $n$ you will probably find it hard to distinguish their effects from an exponential distribution with a parameter close to $\frac{n}{t}$. The distribution is asymptotically exponential because $\left(1-\frac{y}{m}\right)^{m} \to \exp(-y)$ as $m \to \infty$.

Answer 4

To say, "suppose you randomly divide 500 million in income among 10,000 people" is insufficiently specific to answer the question. There are many different random process that could be used to allocate a fixed amount of money to a fixed number of people, and each will have its own characteristics for the resulting distribution. Here are three generative processes I could think of, and the distributions of wealth each creates.

library(MASS)

w <- 500000000 #wealth
p <- 10000 #people

Method 1, posted by OP:

Choose 'p' numbers from [0,w) uniformly at random. Sort these. Append '0' to the front. Hand out dollar amounts represented by the differences between successive elements in this list.

d <- diff(c(0,sort(runif(p-1,max=w)),w)) #wealth-distribution
h <- hist(d, col="red", main="Exponential decline", freq = FALSE, breaks = 45,
     xlim = c(0, quantile(d, 0.99)))
fit <- fitdistr(d,"exponential")
curve(dexp(x, rate = fit$estimate), col = "black", type="p", 
      pch=16, add = TRUE)

Method 2:

Chose 'p' numbers from [0, w) uniformly at random. Consider these 'weights', so 'w' doesn't actually matter at this stage. Normalize the weights. Hand out dollar amounts represented by the fraction of 'w' corresponding to each weight.

d <- runif(p,max=w) #weigh-distribution
d <- d/sum(d)*w #wealth-distribution
h <- hist(d, col="red", main="pretty uniform", freq = FALSE, breaks = 45, 
          xlim = c(0, quantile(d, 0.99)))

Method 3:

Start with 'p' 0s. w times, add 1 to a one of them, selected uniformly at random.

d <- rep(0, p)
for( i in 1:5000000){ ## for-loops in R are terrible, but this gives the idea.
    k <- floor(runif(1, max=p)) + 1    
    d[k] = (d[k] + 1)
}
h <- hist(d, col="red", main="kinda normalish?", freq = FALSE, breaks = 45,
          xlim = c(0, quantile(d, 0.99)))

Answer 5

Let me add something regarding your addendum.

In the continuous case, as pointed out by Glen_b and Henry, the exact PDF for the amount each person receives is \begin{equation} p(x) = \frac{N-1}{X}\left(1-\frac{x}{X}\right)^{N-2}, \end{equation} where $N$ is the number of people, and $X$ is the total amount of money.

In the discrete case, assuming that there are $M$ coins to distribute, the probability for a particular person to receive $m$ coins is \begin{equation} p(m) = \frac{N-1}{M+1}\prod_{j=0}^{N-3}\left(1-\frac{m}{M-j}\right)^{N-2}. \end{equation} When $M\gg N$, two cases agree with each other. For sufficiently large $N$ and as long as we stay away from the tail, they look like exponential distributions.

In both cases, as we are sampling $N$ times from this true probability distribution, there will be error associated with the finite sample size.

However, performing the error analysis does not seem to be straightforward because different samplings in this case are not independent. They have to sum up to the total amount, and how much the first person receives affects the probability distribution for the second person, and so on.

My previous answer does not suffer from this issue, but I think it would be helpful to see how it can be resolved in this approach.

Answer 6

Good theoretical analysis done by the upvoted answers. However, here's my simple, empirical view on why the distribution is exponential.

When you distribute the money randomly, let's consider you do it one-by-one. Let S be the original sum.

For the first man, you must choose a random amount between 0 and S. Thus, on average, you will choose S/2 and remain with S/2.

For the second man, you would choose randomly between 0 and, on average, S/2. Thus, on average, you'll choose S/4 and remain with S/4.

So, you would basically be splitting the sum in half each time (statistically speaking).

Although in a real-life example you will not have continuously halved values, this shows why one should expect the distribution to be exponential.