The “pepperoni pizza problem”

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:
enter image description here

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:enter image description 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))

Answer

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$
enter image description here

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)$:

h(x,t) :=   acos((x^2+1-t^2)/(2*x))+t^2*acos((x^2-1+t^2)/(2*x*t))
  -sqrt(((x+t)^2-1)*(1-(t-x)^2))/2 $

j(n,t) := n * ( (1-1/t)^2*(1-1/t^2)^(n-1) 
  + (2/t^2) *  quad_qag(x * (1-h(x,t)/(%pi*t^2))^(n-1),x,t-1,t ,3)[1]) $

j0(n,t) := n *(1-1/t^2)^(n-1)$


tt : 1/(2*0.1581) $
j(11,tt);
4.521719308511862
j(12,tt);
4.522913706608645
j(13,tt);
4.494540361913981

tt : 1/(2*0.1) $
j(27,tt);
10.37509984083333
j(28,tt);
10.37692859747294
j(29,tt);
10.36601271146961

fpprintprec: 4$
nn : makelist(n, n, 2, round(tt^2*1.4))$
jnn : makelist(j(n,tt),n,nn) $ 
j0nn : makelist(j0(n,tt),n,nn) $

plot2d([[discrete,nn,jnn],[discrete,nn,j0nn]], 
    [legend,"exact","approx"],[style,[linespoints,1,2]],[xlabel,"n"],  [ylabel,"j(n)"],
[title,concat("t=",tt, "  r/R=",1/(2*tt))],[y,2,tt^2*0.5]);

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$)

enter image description here

enter image description here

enter image description here

Attribution
Source : Link , Question Author : Ben S. , Answer Author : leonbloy

Leave a Comment