This problem arose in a different context at work, but I have translated it to pizza.
Suppose you have a circular pizza of radius $R$. Upon this disc, $n$ pepperoni will be distributed completely randomly. All pepperoni have the same radius $r$.
A pepperoni is "free" if it does not overlap any other pepperoni.
You are free to choose $n$.
Suppose you choose a small $n$. The chance that any given pepperoni is free are very large. But $n$ is small so the total number of free pepperoni is small. Suppose you choose a large $n$. The chance that any given pepperoni is free are small. But there are a lot of them.
Clearly, for a given $R$ and $r$, there is some optimal $n$ that maximizes the expected number of free pepperoni. How to find this optimum?
Edit: picking the answer
So it looks like leonbloy's answer given the best approximation in the cases I've looked at:
r/R n* by simulation n_free (sim) (R/2r)^2
0.1581 12 4.5 10
0.1 29 10.4 25
0.01 2550 929.7 2500
(There's only a few hundred trials in the r=0.01 sim, so 2550 might not be super accurate.) So I'm going to pick it for the answer. I'd like to thank everyone for their contributions, this has been a great learning experience.
Here are a few pictures of a simulation for r/R = 0.1581, n=12:
Edit after three answers posted:
I wrote a little simulation. I'll paste the code below so it can be checked (edit: it's been fixed to correctly pick points randomly on a unit disc). I've looked at two three cases so far. First case, r = 0.1581, R = 1, which is roughly p = 0.1 by mzp's notation. At these parameters I got n* = 12 (free pepperoni = 4.52). Arthur's expression did not appear to be maximized here. leonbloy's answer would give 10. I also did r = 0.1, R = 1. I got n* = 29 (free pepperoni = 10.38) in this case. Arthur's expression was not maximized here and leonbloy's answer would give 25. Finally for r = 0.01 I get roughly n*=2400 as shown here:
Here's my (ugly) code, now edited to properly pick random points on a disc:
from __future__ import division
import numpy as np
# the radius of the pizza is fixed at 1
r = 0.1 # the radius of the pepperoni
n_to_try = [1,5,10,20,25,27,28,29,30,31,32,33,35] # the number of pepperoni
trials = 10000# the number of trials (each trial randomly places n pepperoni)
def one_trial():
# place the pepperoni
pepperoni_coords = []
for i in range(n):
theta = np.random.rand()*np.pi*2 # a number between 0 and 2*pi
a = np.random.rand() # a number between 0 and 1
coord_x = np.sqrt(a) * np.cos(theta) # see http://mathworld.wolfram.com/DiskPointPicking.html
coord_y = np.sqrt(a) * np.sin(theta)
pepperoni_coords.append((coord_x, coord_y))
# how many pepperoni are free?
num_free_pepperoni = 0
for i in range(n): # for each pepperoni
pepperoni_coords_copy = pepperoni_coords[:] # copy the list so the orig is not changed
this_pepperoni = pepperoni_coords_copy.pop(i)
coord_x_1 = this_pepperoni[0]
coord_y_1 = this_pepperoni[1]
this_pepperoni_free = True
for pep in pepperoni_coords_copy: # check it against every other pepperoni
coord_x_2 = pep[0]
coord_y_2 = pep[1]
distance = np.sqrt((coord_x_1 - coord_x_2)**2 + (coord_y_1 - coord_y_2)**2)
if distance < 2*r:
this_pepperoni_free = False
break
if this_pepperoni_free:
num_free_pepperoni += 1
return num_free_pepperoni
for n in n_to_try:
results = []
for i in range(trials):
results.append(one_trial())
x = np.average(results)
print "For pizza radius 1, pepperoni radius", r, ", and number of pepperoni", n, ":"
print "Over", trials, "trials, the average number of free pepperoni was", x
print "Arthur's quantity:", x* ((((1-r)/1)**(x-1) - (r/1)) / ((1-r) / 1))
Updated: see below (Update 3)
Here's another approximation. Consider the center of the disks (pepperoni) as an homogeneous point process of density $\lambda = n/A$, where $A=\pi R^2$ is the pizza surface. Let $D$ be nearest neighbour distance from a given center. Then
$$P(D\le d) = 1- \exp(-\lambda \pi d^2)=1- \exp(-n \,d^2/R^2) \tag{1}$$
A pepperoni is free if $D > 2r$. Let $E$ be the expected number of free peperoni.
Then $$E= n\, P(D>2r) = n \exp (-n \,4 \, r^2/R^2) = n \exp(-n \, p)$$ where $p=(2r)^2/R^2$ (same notation as mzp's answer).
The maximum is attained for (ceil or floor of) $n^*=1/p=(R/2r)^2\tag{2}$
Update 1: Formula $(1)$ could be corrected for the border effects, the area near the border would be computed as the intersection of two circles. It looks quite cumbersome, though.
Update 2: In the above, I assumed that the center of the pepperoni could be placed anywhere inside the pizza circle. If that's not the case, if the pepperoni must be fully inside the pizza, then $R$ should be replaced by the "effective radius" $R' = R-r$
Update 3: The Poisson approach is really not necessary. Here's an exact solution
Let $$t = \frac{R}{2r}$$
(Equivalently, think of $t$ as the pizza radius, and assume a pepperoni of radius $1/2$). Assume $t>1$. Let $g(x)$ be the area of a unit circle, at a distance $x$ from the origin, intersected with the circle of radius $t$. Then
$$g(x)=\begin{cases}\pi & {\rm if}& 0\le x \le t-1\\ h(x) & {\rm if}& t-1<x \le t \\ 0 & {\rm elsewhere} \end{cases} \tag{3}$$ where $$h(x)=\cos^{-1}\left(\frac{x^2+1-t^2}{2x}\right)+t^2 \cos^{-1}\left(\frac{x^2+t^2-1}{2xt}\right) -\frac{1}{2}\sqrt{[(x+t)^2-1][1-(t-x)^2]} \tag{4}$$
Here's a graph of $g(x)/\pi$ for $t=5$
Let the random variable $e_i$ be $1$ if the $i$ pepperoni is free, $0$ otherwise. Then
$$E[e_i \mid x] = \left(1- \frac{g(x)}{\pi \, t^2}\right)^{n-1} \tag{5}$$ (remaining pepperoni fall in the free area). And
$$E[e_i] =E[E(e_i \mid x)]= \int_{0}^t \frac{2}{t^2} x \left(1- \frac{g(x)}{\pi \, t^2}\right)^{n-1} dx = \\ =\left(1- \frac{1}{t^2}\right)^{n-1} \left(1- \frac{1}{t}\right)^2 +\frac{2}{t^2} \int_{t-1}^t x \left(1- \frac{h(x)}{\pi\, t^2}\right)^{n-1} dx \tag{6}$$
The objective function (expected number of free pepperoni) is then given by:
$$J(n)=n E[e_i ] \tag{7} $$
This is exact... but (almost?) intractable. However, it can be evaluated numerically [**].
We can also take as approximation $$g(x)=\pi$$ for $0\le x < t$ (neglect border effects) and then it gets simple:
$$E[e_i ] =E[e_i \mid x]= \left(1- \frac{1}{t^2}\right)^{n-1}$$
$$J(n)= n \left(1- \frac{1}{t^2}\right)^{n-1} \tag{8}$$
To maximize, we can write $$\frac{J(n+1)}{J(n)}= \frac{n+1}{n} \left(1- \frac{1}{t^2}\right)=1 $$ which gives
$$ n^{*}= t^2-1 = \frac{R^2}{4 r^2}-1 \tag{9}$$ quite similar to $(2)$.
Notice that as $t \to \infty$, $J(n^{*})/n^{*} \to e^{-1}$, i.e. the proportion of free pepperoni (when using the optimal number) is around $36.7\%$. Also, the "total pepperoni area" is 1/4 of the pizza.
[**] Some Maxima code to evaluate (numerically) the exact solution $(7)$:
The first result agrees with the OP simulation, the second is close: I got $n^{*}=28$ instead of $29$. For the third case, I get $n^{*}=2529$ ($j(n)=929.1865331$)