The Cake Cutting Problem

Introduction

Alon Amit’s post gave a very good explanation about the cake cutting problem. In mathematical language, the problem is stated as:

\begin{align} X_1 & \sim U(0, 1),\\ X_2 & \sim U(0, 1-X_1),\\ X_3 & \sim U(0, 1-X_1-X_2),\\ & \vdots \\ \end{align}

What is the probability that $X_1$ is the maximum among $X_1,\ X_2,\ X_3,\ \ldots$?

The answer is the Golomb-Dickman constant $0.6243\cdots$. Alon Amit made efforts to explain why “the Golomb-Dickman pops up”. In Michael Lamar’s post, there is Matlab/Octave code for a simulation study.

I wrote the following R program:

simu_exp_1 <- function()
{x <- 1
 first_cut <- runif(1, min = 0, max = x)
 temp <- first_cut
 while(1) {
   x <- x - temp
   if(x < first_cut) return(1)
   next_cut <- runif(1, min = 0, max = x)
   if(next_cut > first_cut) return(0)
   temp <- next_cut
 }
}

N <- 1000000
simu_re <- replicate(n = N, simu_exp_1())
mean_re <- cumsum(simu_re) / (1:N)
(mean_re[N])

## [1] 0.623832

Inspired by the above cake cutting problem, I asked myself the following question:

\begin{align} X_1 & \sim U(0, 1),\\ X_2 & \sim U(0, 1-X_1),\\ & \vdots \\ X_n & \sim U(0, 1 - \sum_{i=1}^{n-1} X_i), \end{align}

Define

$$ P_n = \hbox{Prob}\left(X_1 > X_2 > \cdots > X_n\right),\ \hbox{for}\ n = 2,\ 3, \ldots. $$

When $n$ is given, what is $P_n$?

Finding $P_2$

In the case $n=2$, we can find the the joint probability density function of $X_1$ and $X_2$ is

$$ f(x_1, x_2) = I_{(0, 1)}(x_1) I_{(0, 1-x_1)}(x_2)\frac{1}{1-x_1}. $$

Thus,

\begin{align} P_2 &= \int_{x_1> x_2} f(x_1, x_2)dx_1 dx_2\\ &= \int_{0}^1 \int_{0}^{\hbox{min} \{x_1, 1-x_1\}} \frac{1}{1-x_1}dx_2 dx_1\\ &= \int_{0}^1 \frac{1}{1-x_1} \hbox{min} \{x_1, 1-x_1\} dx_1\\ &= \int_{0}^{1/2}\frac{x}{1-x} dx + \int_{1/2}^1 \frac{1-x}{1-x} dx\\ &= \ln(2) \approx 0.6931472 \end{align}

Finding $P_n$ by simulation

My R code is given below

simu_exp_2 <- function(n = 2)
{x <- rep(0, n)
 x[1] <- runif(1, min = 0, max = 1)
 for(i in 2:n) {
   x[i] <- runif(1, min = 0, max = 1 - sum(x[1:(i-1)]))
 }
 y <- sort(x, decreasing = TRUE)
 if(identical(x, y)) return(1)
 return(0)
}

get_result <- function(n = 2) 
{N <- 1000000
 simu_re <- replicate(n = N, simu_exp_2(n))
 mean_re <- cumsum(simu_re) / (1:N)
 mean_re[N]
}

(get_result(n = 2))

## [1] 0.692733

(get_result(n = 3))

## [1] 0.431809

(get_result(n = 4))

## [1] 0.260376