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{array}{rl}{X}_1& \sim U(0,1),\\ X_2& \sim U(0,1-X_1),\\ X_3& \sim U(0,1-X_1-X_2),\\ & ⋮\end{array}$$

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

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{array}{rl}{X}_1& \sim U(0,1),\\ X_2& \sim U(0,1-X_1),\\ & ⋮\\ X_n& \sim U(0,1-\sum _{i=1}^{n-1}{X}_i),\end{array}$$

Define

$$P_n=\text{Prob}(X_1>X_2>\cdots >X_n),\text{for}n=2,3,\dots .$$

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{array}{rl}{P}_2& ={\int }_{x_1>x_2}f(x_1,x_2)d{x}_{1}d{x}_2\\ & ={\int }_0^1{\int }_0^{\text{min}\{x_1,1-x_1\}}\frac{1}{1-x_1}d{x}_{2}d{x}_1\\ & ={\int }_0^1\frac{1}{1-x_1}\text{min}\{x_1,1-x_1\}d{x}_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{array}$$

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