If there are 200 typographical errors randomly distributed in a 500 page manuscript, find the probability that a given page contains exactly 3 errors.

We can abstract this type of problems as follows:

Suppose there are nn distinguishable boxes and kk indistinguishable balls. Now, we randomly put the balls into the boxes. For each of the boxes, what is the probability that it contains mm balls?

For example, if the first page contains 3 errors, the second page contains 197 errors, and the rest of the pages contain no errors, then the situation corresponds to the situation where the first box contains 3 balls, the second box contains 197 balls, and the rest of the boxes contain no balls. The balls are indistinguishable because we can only determine how many errors are on each page but not which errors are on the page.

To deal with the problem, we simply need to find these two numbers:

  • the number of ways to put kk indistinguishable balls into nn distinguishable boxes, and
  • the number of ways to put k−mk-m indistinguishable balls into n−1n-1 distinguishable boxes.

The latter corresponds to the number of ways to put the balls into the boxes provided that we already know that the given box contains mm balls. After we find these two numbers, their ratio is the probability in question.

To find the number of ways to put kk indistinguishable balls into nn distinguishable boxes, we can use the stars and bars method. To see this, we write a special example. Here is an example of n=4n=4 and k=6k=6: ∣⋆⋆∣⋆∣⋆⋆⋆,{}|{}\star{}\star{}|{}\star{}|{}\star{}\star{}\star{}, which corresponds to the distribution 0,2,1,30,2,1,3. We can see that there are n−1n-1 bars and kk stars. Therefore, the number of ways to put the balls is the same as the number of ways to choose the kk positions of the stars among n+k−1n+k-1 positions. Therefore, the number of ways is Nn,k=(n+k−1k)=(n+k−1)!k!(n−1)!.N_{n,k}=\binom{n+k-1}{k}=\frac{\left(n+k-1\right)!}{k!\left(n-1\right)!}. Therefore, the final probability of the given box containing mm balls is Pn,k(m)=Nn−1,k−mNn,k=(n−1)k!(n+k−m−2)!(k−m)!(n+k−1)!.P_{n,k}(m)=\frac{N_{n-1,k-m}}{N_{n,k}} =\frac{\left(n-1\right)k!\left(n+k-m-2\right)!}{\left(k-m\right)!\left(n+k-1\right)!}.

Another easy way to derive this result is by using the generating function. The number Nn,kN_{n,k} is just the coefficient of xkx^k in the expansion of the generating function (1+x+x2+⋯ )n\left(1+x+x^2+\cdots\right)^n. The generating function is just (1−x)−n\left(1-x\right)^{-n}, which can be easily expanded by using the binomial theorem.

We are now interested in the limit n,k→∞n,k\to\infty with λ≔k/n\lambda\coloneqq k/n fixed. By Stirling’s approximation, we have Pn,k(m)∼(n−1)kk+1/2(n+k−m−2)n+k−m−2+1/2(k−m)k−m+1/2(n+k−1)n+k−1+1/2ek−m+n+k−1−k−n−k+m+2.P_{n,k}(m)\sim\left(n-1\right) \frac{k^{k+1/2}\left(n+k-m-2\right)^{n+k-m-2+1/2}}{\left(k-m\right)^{k-m+1/2}\left(n+k-1\right)^{n+k-1+1/2} } \mathrm e^{k-m+n+k-1-k-n-k+m+2}. The 1/21/2’s in the exponents can just be dropped because you may find that if we extract the 1/21/2’s, the factor tends to unity. The exponential is just constant e\mathrm e. Therefore, we have Pn,k(m)∼(n−1)(λn)λn(n+λn−m−2)n+λn−m−2(λn−m)λn−m(n+λn−1)n+λn−1e=(n+λn−m−2n+λn−1)n((n+λn−m−2)λn(λn−m)(n+λn−1))λn(λn−mn+λn−m−2)m(n−1)(n+λn−1)(n+λn−m−2)2e→e−m+1λ+1 em e−m+1λ+1λ(λλ+1)m1λ+1e=(λλ+1)m1λ+1.\begin{align*} P_{n,k}(m)&\sim\left(n-1\right) \frac{\left(\lambda n\right)^{\lambda n}\left(n+\lambda n-m-2\right)^{n+\lambda n-m-2} } {\left(\lambda n-m\right)^{\lambda n-m}\left(n+\lambda n-1\right)^{n+\lambda n-1}}\mathrm e\\ &=\left(\tfrac{n+\lambda n-m-2}{n+\lambda n-1}\right)^n \left(\tfrac{\left(n+\lambda n-m-2\right)\lambda n}{\left(\lambda n-m\right)\left(n+\lambda n-1\right)}\right)^{\lambda n} \left(\tfrac{\lambda n-m}{n+\lambda n-m-2}\right)^m \tfrac{\left(n-1\right)\left(n+\lambda n-1\right)}{\left(n+\lambda n-m-2\right)^2}\mathrm e\\ &\to\mathrm e^{-\frac{m+1}{\lambda+1}}\,\mathrm e^m\, \mathrm e^{-\frac{m+1}{\lambda+1}\lambda}\left(\tfrac\lambda{\lambda+1}\right)^m\tfrac1{\lambda+1}\mathrm e\\ &=\left(\tfrac\lambda{\lambda+1}\right)^m\tfrac1{\lambda+1}. \end{align*} This is just the geometric distribution with parameter p=1/(λ+1)=n/(k+n)p=1/(\lambda+1)=n/(k+n).

If you want to simulate the number of balls in a box, here is a simple way to do this. First, because each box is the same, we can just focus on the first box without loss of generality. Then, we just need to randomly generate the positions of the n−1n-1 bars among the n+k−1n+k-1 positions, and then return the index of the first bar (which is the number of balls in the first box).

We can then write the following Ruby code to simulate the number of balls in the first box:

def simulate n, k
  (n-1).times.inject(npkm1 = n+k-1) { |bar, i| [rand(npkm1 - i), bar].min }
end

Compare the simulated result with the theoretical result:

def frequency m, n, k, trials
  trials.times.count { simulate(n, k) == m } / trials.to_f
end

def truth m, n, k
  (n-1) * (k-m+1..k).reduce(1,:*) / (n+k-m-1..n+k-1).reduce(1,:*).to_f
end

def approx m, n, k
  n*k**m / ((n+k)**(m+1)).to_f
end

srand 1108
m, n, k = 3, 5000, 8000
p frequency m, n, k, 10000 # => 0.0902
p truth m, n, k # => 0.08965012972626446
p approx m, n, k # => 0.08963271594131858