Computing expected moments of the Rényi parking problem on the circle

M. Hegland C.J. Burden Z. Stachurski

Abstract

A highly accurate and efficient method to compute the expected values of the count, sum, and squared norm of the sum of the centre vectors of a random maximal sized collection of non-overlapping unit diameter disks touching a fixed unit-diameter disk is presented. This extends earlier work on Rényi’s parking problem [Magyar Tud. Akad. Mat. Kutató Int. Közl. 3 (1-2), 1958, pp. 109-127]. Underlying the method is a splitting of the the problem conditional on the value of the first disk. This splitting is proven and then used to derive integral equations for the expectations. These equations take a lower block triangular form. They are solved using substitution and approximation of the integrals to very high accuracy using a polynomial approximation within the blocks.

1 Introduction

Consider the two-dimensional process of unit-diameter disks attaching randomly onto the circumference of a fixed central unit-diameter disk, as shown in Figure 1(a). The accretion is assumed to occur sequentially in such a way that the location of each additional disk is randomly and uniformly distributed over the available accessible part of the central disk. The process stops when no further attachment points are available. Herein we analyse the distributional properties of the vector sum of the location of the centres of the attached disks.

The problem is related to a well-known one-dimensional packing problem known as Rényi’s car parking problem [2]. [8] (see [5, pp 203–218] for an English translation) analysed a stochastic process arising when sequentially packing unit intervals $[Y_{i},Y_{i}+1)$ representing parked cars into an interval $[a,b]$ (where $a,b\in\mathbb{R}$ and $Y_{i}\in[a,b-1]$ is a random number), such that the unit intervals do not overlap. The location of each additional car is randomly and uniformly distributed over the available subset of $[a,b]$ , and the packing is completed when all the remaining gaps are less than one. Rényi showed [8] that the expected value of the number of parked cars satisfies an integral equation resulting from a splitting property of the number of cars in sub-intervals on either side of the first parked car. He further established an exact value (approximately equal to $0.75$ ) for the limit of the expected packing density as the parking interval $b-a\to\infty$ . After the first external disk has attached, our approach involves of mapping the available part of the circumference of the central disk onto the interval $[1,6]$ , as shown in Figure 1 (b) for the interval case which is appropriately mapped to the circular case in Figure 1 (a), and developing analogues of Rényi’s integral equation for features relevant to the disk accretion problem.

The expected maximal number of parked cars is computed explicitly for small $x=b-a$ . However the formulas get very complicated even for moderate $x$ , take a long time to compute and often evaluate poorly due to rounding errors. Thus numerical approaches were established for computing the expectation [3, 4]. These numerical approaches give highly accurate results. The approach taken here is similar in that it also uses a highly accurate piecewise polynomial approximation. In addition to computing the expected number of added disks, we compute the mean squared shift of the vector sum of centres of the attached disks. This requires the solution to high numerical accuracy of a system of integral equations. We review the underlying splitting property in Section 2 for various features relevant to the disk accretion problem and derive a system of integral equations for their expected values in Section 3. The method is then implemented using the new framework developed by [7]. The algorithm is summarised in Section 4 and final results are presented in Section 5.

Refer to caption — Figure 1: (a) The accretion of disks to a central disk, and (b) the equivalent Rényi parking problem, in which the centre of each disk is mapped to the left hand end of a unit-length car. The vector-valued function $v$ is defined by Equation (3). Without loss of generality, the first attached disk is assumed to be centred at $v(0)=(1,0)^{\rm T}$ . Equivalently, the first car is assumed to be located at position 0, and hence the Rényi process pertains to the interval $[1,6]$ .

2 Mathematical modelling

2.1 Rényi Point Process for $[a,b]$

Let $K(a,b)$ be the (random) number of unit-length cars parked in an interval $[a,b]$ when packing is completed. Rényi’s stochastic process defines a random set

\alpha(a,b)=\{Y_{1},\ldots,Y_{K(a,b)}\}\subset[a,b-1],

which is a translational invariant point process, as shown below where the random variables $Y_{i}$ are formally introduced.

Lemma 1 (Translational Invariance).

\alpha(a+t,b+t)=t+\alpha(a,b),\;\text{for all $t\in\mathbb{R}$}.

Proof idea:

The stochastic process for the interval $[a+t,b+t]$ generates a sequence of random variables $Y_{1}+t,\ldots,Y_{K}+t$ where the $Y_{i}$ are the random variables generated by the process for the interval $[a,b]$ . $\spadesuit$ It follows that $K(a,a+x)=K(0,x)$ for all real $x\geq 0$ .

A consequence of the one-dimensionality of the Rényi point process is a splitting property of the random set conditional on the first point $Y_{1}$ . This conditional set is denoted by $\alpha(a,b\mid Y_{1}=y)$ . The splitting property gives rise to the integral equations derived in the next section.

Lemma 2 (Splitting Property).

\alpha(a,b\mid Y_{1}=y)=\{y\}\cup\alpha(a,y)\cup\alpha(y+1,b).

Proof:

When the random variables $Y_{2},Y_{3},\ldots$ are generated, each element is either in $[a,y-1]$ or in $[y+1,b-1]$ . Now $Y_{i},Y_{j}$ are independent if $Y_{i}\in[y+1,b-1]$ and $Y_{j}\in[a,y-1]$ as they cannot overlap. Consequently, the points in the two subsets generated for the conditional process produce two independent random sets $\alpha(a,y)$ and $\alpha(y+1,b)$ . $\spadesuit$

The Rényi point process gives rise to a random counting measure with parameters $a,b$ defined by

N(a,b)=\sum_{y\in\alpha(a,b)}\delta_{y}

where $\delta_{y}$ is the atomic measure for which

\delta_{y}(M)=\begin{cases}1,&\text{if $y\in M$}\\ 0,&\text{otherwise}\end{cases}

for any set $M\subset\mathbb{R}$ . It follows that the support of this measure is $\alpha(a,b)$ . This measure is defined for any set $M\subset\mathbb{R}$ by

N(a,b)(M)=\lvert\alpha(a,b)\cap M\rvert.

(1)

In particular, one has $K(a,b)=N(a,b)([a,b])$ . Note that $K$ is translational invariant, that is $K(a+t,b+t)=K(a,b)$ .

The conditional counting measure $N(a,b\mid Y_{1}=y)$ satisfies a splitting property

N(a,b\mid Y_{1}=y)=\delta_{y}+N(a,y)+N(y+1,b).

(2)

This is a consequence of Lemma 2.1.

2.2 Randomly packing disks touching a central disk

We map the interval $[0,6)$ bijectively onto the unit circle $S^{2}$ by

v(y)=Q(y)\begin{bmatrix}1\\ 0\end{bmatrix},\quad y\in[0,6).

(3)

where

Q(y)=\begin{bmatrix}\cos(\pi y/3)&-\sin(\pi y/3)\\ \sin(\pi y/3)&\cos(\pi y/3)\end{bmatrix},

(4)

as shown in Figure 1. Packing subintervals of $[1,6)$ then corresponds to packing unit diameter disks touching a central unit diameter disk, conditional on an initial disk being attached with its centre at $v(0)=(1,0)^{\rm T}$ . The centres of the touching disks are on the unit-radius circle $S^{2}$ . One can thus apply Rényi’s point processes to the addition of disks subsequent to the first added disk.

We now introduce real and vector valued features of the point process which are polynomial functions of the random variables $v(Y_{1}),\ldots,v(Y_{K})$ invariant under permutations of the $Y_{i}$ . In particular we consider features invariant under any translation of the parameters $a$ and $b$ . Clearly, $K(a,b)$ is an example as $K(a+t,b+t)=K(a,b)$ and $K(a,b)$ is a constant function of the $Y_{i}$ . We introduce a 2-component vector-valued polynomial

F(a,b)=\sum_{i=1}^{K(a,b)}v(Y_{i}).

(5)

The value $((1,0)^{\rm T}+F(1,6))/(1+K(1,6))$ is the centre of mass of the circular polygon with vertices $v(Y_{i})$ . The splitting property of $F(a,b)$ follows from the definition and is

F(a,b\mid Y_{1}=y)=F(a,y)+F(y+1,b)+v(y).

(6)

But $F$ is not invariant under translations since

F(a+t,b+t)=Q(t)F(a,b).

(7)

We introduce the feature

L^{2}(a,b)=\lVert F(a,b)\rVert^{2}

(8)

which is used in Section 5 to determine the expected shift in the centre of the polygon with vertices $v(Y_{i})$ . We introduce another second degree feature

E_{2}(a,b)=\sum_{Y_{i}<Y_{j}}v(Y_{i})^{T}v(Y_{j})=\sum_{Y_{i}<Y_{j}}\cos(\pi(Y_{i}-Y_{j})/3)

(9)

which is better suited to the computations performed in Section 3. We show below how this feature relates to $L^{2}(a,b)$ . This feature is also invariant under translations of the $Y_{i}$ as it only depends on the differences $Y_{i}-Y_{j}$ and the translated difference is $(Y_{i}+t)-(Y_{j}+t)=Y_{i}-Y_{j}$ . This feature admits the following splitting property:

Proposition 3.

\begin{split}E_{2}(a,b\mid Y_{1}=y)=E_{2}(a,y)+E_{2}(y+1,b)+F(a,y)^{T}F(y+1,b)\\ +v(y)^{T}(F(a,y)+F(y+1,b))\end{split}

(10)

Proof:

The sum in Equation (9) decomposes into three sums over $i,j\neq 1$ , one sum over $i$ where $j=1$ and one sum over $j$ where $i=1$ as follows

\sum_{Y_{i}<Y_{j}}=\sum_{Y_{i}<Y_{j}<y}+\sum_{y<Y_{i}<Y_{j}}+\sum_{Y_{i}<y<Y_{j}}+\sum_{Y_{1}=y<Y_{j}}+\sum_{Y_{i}<y=Y_{1}}.

The summand in each term is $\cos(\pi(Y_{i}-Y_{j})/3)$ . $\spadesuit$ We use the splitting property to derive integral equations in the next section for the expectation of $F(a,b)$ and of $E_{2}(a,b)$ .

We now have three features $K,L^{2}$ and $E_{2}$ connected by the following Lemma.

Proposition 4.

L^{2}(a,b)=K(a,b)+2E_{2}(a,b).

Proof:

One has

	$\displaystyle L^{2}(a,b)$	$\displaystyle=\sum_{i=1}^{K(a,b)}\sum_{j=1}^{K(a,b)}v(Y_{i})^{T}v(Y_{j})$
		$\displaystyle=\sum_{i=1}^{K(a,b)}\lVert v(Y_{i})\rVert^{2}+2\sum_{Y_{i}<Y_{j}}v(Y_{i})^{T}v(Y_{j}).$

The result follows as $v(y)$ lies on the unit circle. $\spadesuit$

3 Integral equations

3.1 Equation for $\mathbb{E}K(0,x)$

The random variables $K(a,b)$ denoting the counts of parked vehicles (or disks) in Rényi’s model on an interval $[a,b]$ are translational invariant, that is $K(a,b)=K(a+t,b+t)$ for any $t\in\mathbb{R}$ . An application of these properties provides the following lemma.

Proposition 5.

Let $u_{1}(x)=\mathbb{E}K(0,x)$ be the expected number of vehicles in the interval $[0,x]$ . Then

u_{1}(x)=\frac{2}{x-1}\int_{0}^{x-1}u_{1}(y)\,dy+1.

(11)

Proof:

Let $Y_{1}$ be the (random) position of the first vehicle. For any $a\leq y<b$ , Equation (2) implies that the counting variable $K(a,b)$ conditional on $Y_{1}=y$ satisfies a splitting property

K(a,b\mid Y_{1}=y)=K(a,y)+K(y+1,b)+1.

As $Y_{1}$ is uniformly distributed over $[a,b-1]$ one has

\mathbb{E}K(a,b)=\frac{1}{b-1-a}\int_{a}^{b-1}\mathbb{E}K(a,b\mid Y_{1}=y)\,dy.

As the counting variable is translational invariant, that is $K(a,b)=K(0,b-a)$ one solves the proposition by substituting the conditional count on the right-hand side using the splitting property and substituting $a=0$ and $b=x$ . $\spadesuit$

The expected count for a general interval is

\mathbb{E}K(a,b)=u_{1}(b-a).

3.2 Equation for $\mathbb{E}F(0,x)$

Proposition 6.

Let the (vector valued) function $u_{2}(x)=\mathbb{E}F(0,x)$ where $F(0,x)$ is defined in Equation (5). Then

u_{2}(x)=\frac{1}{x-1}\int_{0}^{x-1}\left(u_{2}(y)+Q(x-y)u_{2}(y)+v(y)\right)\,dy

(12)

where $Q$ is defined in Equation (4).

Proof:

The splitting condition (6) is

F(0,x\mid Y_{1}=y)=F(0,y)+F(y+1,x)+v(y).

Using the translation formula Equation (7) for $F$ and taking the expectation gives

u_{2}(x\mid Y_{1}=y)=u_{2}(y)+Q(y+1)u_{2}(x-1-y)+v(y).

Integrating over $y$ and reversing the order of the integration for the second term gives the result. $\spadesuit$

3.3 Equation for $\mathbb{E}E_{2}(0,x)$

Proposition 7.

Let $u_{3}(x)=\mathbb{E}E_{2}(0,x)$ . Then

u_{3}(x)=\frac{2}{x-1}\int_{0}^{x-1}u_{3}(y)\,dy+g_{3}(x)

(13)

where

\begin{split}g_{3}(x)=\frac{1}{x-1}\int_{0}^{x-1}\left(u_{2}(y)^{T}Q(y+1)u_{2}(x-1-y)+v(y)^{T}u_{2}(y)\right.\\ \left.+\,v(-1)^{T}u_{2}(y)\right)\,dy\end{split}

Proof:

The splitting of $E_{2}$ , Equation (10) gives

\begin{split}E_{2}(0,x\mid Y_{1}=y)=E_{2}(0,y)+E_{2}(y+1,x)+F(0,y)^{T}F(y+1,x)\\ +\,v(y)^{T}\left(F(0,y)+F(y+1,x)\right).\end{split}

Using the independence of $F(0,y)$ and $F(y+1,x)$ one gets

\begin{split}u_{3}(x\mid Y_{1}=y)=u_{3}(y)+u_{3}(x-1-y)+u_{2}(y)^{T}Q(y+1)u_{2}(x-1-y)\\ +\,v(y)^{T}\left(u_{2}(y)+Q(y+1)u_{2}(x-1-y)\right).\end{split}

Noting that $v(y)^{T}Q(y+1)=v(-1)^{T}$ and integrating over $y$ with a change of integration variable in the last term gives the result $\spadesuit$

4 Numerical solution

Equations (11), (12) and (13), for the expectations $u_{1}(x)=\mathbb{E}K(0,x)$ , $u_{2}(x)=\mathbb{E}F(0,x)$ and $u_{3}(x)=\mathbb{E}E_{2}(0,x)$ are all integral equations of the form

u=\mathcal{L}u+g

(14)

for some linear integral operator $\mathcal{L}$ and function $g$ specified in the following sub-sections. In the Rényi process, a space of length less than 1 must be unoccupied, so $u(x)=0$ for $x\in[0,1)$ . Furthermore,

\mathcal{L}u(x)=\frac{1}{x-1}\int_{0}^{x-1}M(x,y)u(y)\,dy,\quad\text{for $x>1$},

where $M(x,y)=I+Q(x-y)$ for the vector case $u_{2}(x)$ , and $M(x,y)=2$ for the scalar cases $u_{1}(x)$ and $u_{3}(x)$ . We are interested in computing $u(x)$ for $x\in[0,5]$ . The solution of Equation (14) follows trivially for $x\in[0,2)$ :

u(x)=\begin{cases}0,&x\in[0,1);\\ g(x),&x\in[1,2).\end{cases}

Now if one interprets the function $u(x)$ as a block vector where each block corresponds to a function on an interval $[k-1,k]$ then the integral equation (14) is of block lower triangular structure and can be solved by repeated substitution. Specifically, one computes $u(x)$ for $x\in[k,k+1]$ by

u(x)=\frac{1}{x-1}\int_{0}^{x-1}M(x,y)u(y)+g(x),\quad x\in[k,k+1)

which uses $u(x)$ for $x\in[k-1,k]$ . It then follows directly that $u$ is $C^{\infty}$ on each interval $[j-1,j]$ for $j=1,2,\ldots$ and is continuous on $[1,\infty)$ .

The approach was implemented in the Julia language [1] using the ApproxFun package [7, 6]. The code is available on request from the first author. Plots of the functions $u_{1}$ , $u_{2}$ and $u_{3}$ are shown in Figure 2. Computed values agree to machine accuracy with analytic calculations of $u_{1}(x)$ for $0\leq x\leq 4$ [9], and our own analytic calculations of $u_{1}(5)$ and $u_{2}(x)$ and $u_{3}(x)$ for $0\leq x\leq 3$ .

5 Results

It remains to compute the expected values of features pertaining to the full set of attached disks. The expected total number of disks is

\mathbb{E}[1+K(1,6)]=1+u_{1}(5)=4.48508592498075.

The expected vector sum of centres of the attached discs conditional on the first disk being located at $(1,0)^{\rm T}$ is

$\displaystyle\mathbb{E}\left[\begin{pmatrix}1\\ 0\end{pmatrix}+F(1,6)\right]$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}1\\ 0\end{pmatrix}+Q(1)\mathbb{E}[F(0,5)]$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}1\\ 0\end{pmatrix}+\begin{pmatrix}\cos(\pi/3)&-\sin(\pi/3)\\ \sin(\pi/3)&\cos(\pi/3)\end{pmatrix}u_{2}(5)$
	$\displaystyle=$	$\displaystyle\begin{pmatrix}\text{0.00226785060421}\\ 0\end{pmatrix}.$

By symmetry, the second component must be zero, and this is confirmed to machine accuracy. Finally, define the mean square shift of the sum of centres of the attached disks as

$\displaystyle I$	$\displaystyle:=$	$\displaystyle\left\lVert\begin{pmatrix}1\\ 0\end{pmatrix}+F(1,6)\right\rVert^{2}$
	$\displaystyle=$	$\displaystyle 1+2(1,0)Q(1)F(0,5)+\lVert F(0,5)\rVert^{2}$
	$\displaystyle=$	$\displaystyle 1+2\left(\cos(\pi/3),-\sin(\pi/3)\right)F(0,5)+L^{2}(0,5).$

Its expected value is, using Proposition 2.2,

\mathbb{E}[I]=1+2\left(\cos(\pi/3),-\sin(\pi/3)\right)u_{2}(5)+u_{1}(5)+2u_{3}(5)=0.2325047203936.

Acknowledgements

An initial version of the code was developed and tested using the cloud service for VSCode on juliahub.com.

References

[1] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B. Shah “Julia: A Fresh Approach to Numerical Computing” In SIAM Review 59.1, 2017, pp. 65–98 DOI: 10.1137/141000671
[2] Matthew P. Clay and Nándor J. Simányi “Rényi’s parking problem revisited” In Stoch. Dyn. 16.2, 2016, pp. 1660006, 11 DOI: 10.1142/S0219493716600066
[3] M. Lal and P. Gillard “Evaluation of a constant associated with a parking problem” In Math. Comp. 28, 1974, pp. 561–564 DOI: 10.2307/2005927
[4] George Marsaglia, Arif Zaman and John C. W. Marsaglia “Numerical solution of some classical differential-difference equations” In Math. Comp. 53.187, 1989, pp. 191–201 DOI: 10.2307/2008355
[5] Institute Mathematical Statistics “Selected Translations in Mathematical Statistics and Probability.” American Mathematical Society, Providence, RI, 1963
[6] S. Olver “ApproxFun.jl v0.13.14”, 2022 URL: https://github.com/JuliaApproximation/ApproxFun.jl
[7] S. Olver and A Townsend “A practical framework for infinite-dimensional linear algebra” In Proceedings of the 1st First Workshop for High Performance Technical Computing in Dynamic Languages, 2014, pp. 57–62
[8] Alfréd Rényi “On a one-dimensional problem concerning random space filling” In Magyar Tud. Akad. Mat. Kutató Int. Közl. 3.1-2, 1958, pp. 109–127
[9] Howard Weiner “Elementary Treatment of the Parking Problem” In Sankhyā: The Indian Journal of Statistics, Series A (1961-2002) 31.4 Springer, 1969, pp. 483–486 URL: http://www.jstor.org/stable/25049616

Computing expected moments of the Rényi parking problem on the circle

Abstract

1 Introduction

2 Mathematical modelling

2.1 Rényi Point Process for [a,b][a,b]

Lemma 1 (Translational Invariance).

Proof idea:

Lemma 2 (Splitting Property).

Proof:

2.2 Randomly packing disks touching a central disk

Proposition 3.

Proof:

Proposition 4.

Proof:

3 Integral equations

3.1 Equation for 𝔼​K​(0,x)\mathbb{E}K(0,x)

Proposition 5.

Proof:

3.2 Equation for 𝔼​F​(0,x)\mathbb{E}F(0,x)

Proposition 6.

Proof:

3.3 Equation for 𝔼​E2​(0,x)\mathbb{E}E_{2}(0,x)

Proposition 7.

Proof:

4 Numerical solution

5 Results

Acknowledgements

References

Author address

2.1 Rényi Point Process for $[a,b]$

3.1 Equation for $\mathbb{E}K(0,x)$

3.2 Equation for $\mathbb{E}F(0,x)$

3.3 Equation for $\mathbb{E}E_{2}(0,x)$