Mean perimeter and mean area of the convex hull over planar random walks

TL;DR

This paper analyzes the geometric properties of convex hulls formed by planar random walks, deriving asymptotic behaviors for mean perimeter and area, with applications in biology and ecology.

Contribution

It provides the first derivation of the large $n$ asymptotics of the mean perimeter and area of convex hulls for planar random walks with symmetric jumps, including finite-size corrections.

Findings

Mean perimeter asymptotics derived for large n

Mean area computed for isotropic Gaussian jumps

Finite-size corrections improve accuracy for small n

Abstract

We investigate the geometric properties of the convex hull over $n$ successive positions of a planar random walk, with a symmetric continuous jump distribution. We derive the large $n$ asymptotic behavior of the mean perimeter. In addition, we compute the mean area for the particular case of isotropic Gaussian jumps. While the leading terms of these asymptotics are universal, the subleading (correction) terms depend on finer details of the jump distribution and describe a "finite size effect" of discrete-time jump processes, allowing one to accurately compute the mean perimeter and the mean area even for small $n$, as verified by Monte Carlo simulations. This is particularly valuable for applications dealing with discrete-time jumps processes and ranging from the statistical analysis of single-particle tracking experiments in microbiology to home range estimations in ecology.

Figures

Equations

$$⟨S_2⟩=\sqrt{16\pi Dt},$$

$$⟨V_2⟩=\pi Dt.$$

$$⟨L_n⟩=2\sum _{k=1}^n\frac{⟨|x_k+i{y}_k|⟩}{k}$$

$$(x_0,y_0)=(0,0),(x_k,y_k)=(x_{k-1},y_{k-1})+({\eta }_k^x,{\eta }_k^y)(k=1,2,\dots ,n),$$

$$\P \{{\eta }_k^x\in (x,x+dx),{\eta }_k^y\in (y,y+dy)\}=p(x,y)dxdy.$$

$$p(-x,-y)=p(x,y).$$

$${\hat{\rho }}_{\theta }(k)=\underset{-\infty }{\overset{\infty }{\int }}dx\underset{-\infty }{\overset{\infty }{\int }}dyp(x,y)e^{ik(x\cos \theta +y\sin \theta )},$$

$${\hat{\rho }}_{\theta }(k)\simeq 1-|a_{\theta }k{|}^{\mu }+o(|k{|}^{\mu })(k\to 0),$$

$$⟨L_n⟩\simeq C_0{n}^{\frac{1}{2}}+C_1+o(1)(n\gg 1),$$

$$C_0=\frac{\sqrt{2}}{\sqrt{\pi }}\underset{0}{\overset{2\pi }{\int }}d\theta {\sigma }_{\theta },C_1=\underset{0}{\overset{2\pi }{\int }}d\theta {\sigma }_{\theta }{\gamma }_{\theta },$$

$${\sigma }_{\theta }^2$$

$${\gamma }_{\theta }$$

$$⟨L_n⟩\simeq C_0{n}^{\frac{1}{2}}+C_1+C_2{n}^{-\frac{1}{2}}+o(n^{-\frac{1}{2}})(n\gg 1),$$

$$C_2=\frac{C_0}{8}+\frac{\sqrt{2}}{24\sqrt{\pi }}\underset{0}{\overset{2\pi }{\int }}d\theta {\sigma }_{\theta }{\mathcal{K}}_{\theta }$$

$${\mathcal{K}}_{\theta }=\frac{1}{{\sigma }_{\theta }^4}\underset{k\to 0}{\lim }\frac{{\partial }^4{\hat{\rho }}_{\theta }(k)}{\partial k^4}=\frac{⟨({\eta }^x\cos \theta +{\eta }^y\sin \theta {)}^4⟩}{⟨({\eta }^x\cos \theta +{\eta }^y\sin \theta {)}^2{⟩}^2}.$$

$${\hat{\rho }}_{\theta }(k)\simeq 1-|a_{\theta }k{|}^{\mu }+b_{\theta }|k{|}^{\nu }+o(|k{|}^{\nu })(k\to 0),$$

$$⟨L_n⟩\simeq C_0{n}^{1/\mu }+C_1{n}^{1-(\nu -1)/\mu }+o(n^{1-(\nu -1)/\mu })(n\gg 1),$$

$$C_0$$

$$C_1$$

$$⟨L_n⟩\simeq C_0{n}^{1/\mu }+C_1+o(1)(n\gg 1),$$

$$C_1=\underset{0}{\overset{2\pi }{\int }}d\theta {\gamma }_{\theta },$$

$${\gamma }_{\theta }=\frac{1}{\pi }\underset{0}{\overset{\infty }{\int }}\frac{dk}{k^2}\ln \left(\frac{1-{\hat{\rho }}_{\theta }(k)}{(ak{)}^{\mu }}\right).$$

$$\gamma =a\frac{\zeta (1/\mu )}{(2\pi {)}^{1/\mu }\sin (\pi /(2\mu ))},$$

$${\sigma }^{-2}⟨A_n⟩=\frac{\pi }{2}n+\gamma \sqrt{8\pi }{n}^{\frac{1}{2}}+\pi (\mathcal{K}/12+{\gamma }^2)+o(1),$$

$$\gamma =\frac{\zeta (1/2)}{\sqrt{2\pi }}=-0.58259\dots$$

$${\sigma }^{-2}⟨A_n⟩=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^{n-i}\frac{1}{\sqrt{ij}}.$$

$$M(\theta)=\max\limits_{0\leq s\leq 1}\bigl{\{}X(s)\cos\theta+Y(s)\sin\theta\bigr{\}}.$$

$$L$$

$$A$$

$$M(\theta)=\max\limits_{0\leq k\leq n}\bigl{\{}x_{k}\cos\theta+y_{k}\sin\theta\bigr{\}}.$$

$$M^{\prime }(\theta )=-x_{k^{\ast }}\sin \theta +y_{k^{\ast }}\cos \theta .$$

$$⟨L_n⟩$$

Full text

Keywords: convex hull, random walk, maximum statistics, diffusion, Lévy flights

1 Introduction

Consider a set of $n$ points with position vectors $\{\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n}\}$ in a $d$-dimensional space. The most natural and perhaps the simplest way to characterize the shape of this set of points is by drawing the convex hull around this set: a convex hull is the unique minimal convex polytope that encloses all the points. This unique polytope is convex since the line segment joining any two points on the surface of the polytope is fully contained within the polytope. Properties of such convex polytopes have been widely studied in mathematics, computer science (image processing and patter recognition) and in the physics of crystallography (the Wulff construction). In two dimension, where the convex hull is a polygon, there are many other applications, most notably in ecology where the home range of animals or the spread of an epidemics are typically estimated by convex hulls. For a review on the history and applications of convex hulls, see Ref. [1].

When the points $\{\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n}\}$ are drawn randomly from a joint distribution $P(\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n})$, the associated convex hull also becomes random and characterizing its statistical properties is a challenging problem, since the convex hull is a highly nontrivial functional of the random variables $\{\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n}\}$. For instance, what can one say about the statistics of the surface area $S_{d}$ or the volume $V_{d}$ of the convex hull, for a given joint distribution $P(\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n})$? Even finding the mean surface area $\langle S_{d}\rangle$ or the mean volume $\langle V_{d}\rangle$, for arbitrary joint distribution, is a formidably difficult problem. In the special case when the points are independent and identically distributed, i.e., when the joint distribution factorizes as $P(\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n})=\prod_{k=1}^{n}P(\vec{r}_{k})$ (with $P(\vec{r}_{k})$ representing the marginal distribution), several results on the statistics of the surface and volume of the convex hull are known (see [1] for a historical review). However, for correlated points where the joint distribution does not factorize, very few results are available.

The simplest example of a set of correlated points corresponds to the case of a random walk in $d$-dimensional continuous space, where $\vec{r}_{k}$ represents the position of the walker at step $k$, starting at the origin at step [math]. The position evolves via the Markov rule, $\vec{r}_{k}=\vec{r}_{k-1}+\vec{\eta}_{k}$, where $\vec{\eta}_{k}$ represents the jump at step $k$, and one assumes that $\vec{\eta}_{k}$ are independent and identically distributed random variables, each drawn from some prescribed distribution $p(\vec{\eta}_{k})$. The walk evolves up to $n$ steps generating the vertices $\{\vec{r}_{1},\vec{r}_{2},\ldots,\vec{r}_{n}\}$ of its trajectory. There is a unique convex hull for each sample of this trajectory and what can one say about the mean surface area or the mean volume of this convex hull, given the jump distribution $p(\vec{\eta}_{k})$? This is the basic problem that we address in this paper. We show that at least for $d=2$ (planar random walks), it is possible to obtain precise explicit results for all $n$ for the mean perimeter and the mean area of the convex hull of the walk, for a large class of jump distributions $p(\vec{\eta})$, including in particular Lévy flights where the jump distribution has a fat tail. We also obtain similar results for the mean area of the convex hull but under additional assumptions on the jump distribution.

This problem concerning the convex hull of a random walk becomes somewhat simpler in the special case of the Brownian limit, where several results are known. Consider, for example a jump distribution $p(\vec{\eta}_{k})$ with zero mean and a finite variance $\sigma^{2}$. In this case, the walk converges in the large $n$ limit to the Brownian motion. In other words, one can consider the continuous-time limit, as $\sigma^{2}\to 0$ and $n\to\infty$ with $n\sigma^{2}=2\,D\,t$ being fixed (here $D$ is called the diffusion constant and $t$ is the duration of the walk). In this Brownian limit and for $d=2$, the mean perimeter and the mean area have been known exactly for a while. Takács [2] computed the mean perimeter

[TABLE]

while El Bachir [3] and Letac [4] computed the mean area

[TABLE]

For a planar Brownian bridge of duration $t$ (where the walker returns to the origin after time $t$), the mean perimeter $\langle S_{2}\rangle_{\rm bridge}=\sqrt{\pi^{3}\,D\,t}$ was computed by Goldman [5], while the mean area $\langle V_{2}\rangle_{\rm bridge}=(2\pi/3)\,D\,t$ was computed relatively recently by Randon-Furling et al. [6]. An interesting extension of this problem in $d=2$ is to the case of $N$ independent planar Brownian motions (or Brownian bridges) [6, 1]. This is relevant in the context of the home range of animals, where $N$ represents the size of an animal population and the trajectory of each animal is approximated by a Brownian motion during their foraging period. For a fixed population size $N$, the mean perimeter and the mean area of the convex hull was computed exactly: $\langle S_{2}\rangle=\alpha_{N}\,\sqrt{D\,t}$ and $\langle V_{2}\rangle=\beta_{N}\,D\,t$, where the prefactors $\alpha_{N}$ and $\beta_{N}$ were found to have nontrivial $N$ dependence [6, 1]. For $d>2$, very few exact results are known for this problem. For a single $(N=1)$ Brownian motion, the mean surface area and the mean volume of the convex hull was recently computed by Eldan [7]: $\langle S_{d}\rangle=\frac{2(4\pi Dt)^{(d-1)/2}}{\Gamma(d)}$ and $\langle V_{d}\rangle=\frac{(\pi Dt)^{d/2}}{\Gamma(d/2+1)^{2}}$ (see also [8] for another derivation and extension to Brownian bridges). However, for $N>1$ and $d>2$, no exact result is available. Finally, going beyond the mean surface and the mean volume, very few results are known for higher moments (see the review [1] for results on variance) or even the full distribution of the surface or the volume of the convex hull of Brownian motion (see Refs. [9, 10] for a recent discussion on the distribution of the perimeter in $d=2$ and $N=1$). Very recently, the full distribution (including the large deviation tails) of the perimeter and the area of $N\geq 1$ planar Brownian motions were calculated numerically [11, 12]. Some rigorous results on the convex hulls of Lévy processes were recently derived [13, 14].

If one is interested only in the mean area or the mean volume of the convex hull of a generic stochastic process (not necessarily just a random walk), a particular simplification occurs in $d=2$ (planar case) where several analytical results can be derived by adapting Cauchy’s formula [15, 16] for arbitrary closed convex curves in $d=2$. Indeed by employing Cauchy’s formula for every realization of a random planar convex hull, it was shown in Refs. [6, 1] that the problem of computing the mean perimeter and the mean area of an arbitrary two dimensional stochastic process (can in general be non-Markovian and in discrete-time) can be mapped to computing the extremal statistics associated with the one dimensional component of the process (see Section 3 for the precise mapping). This mapping was introduced originally in [6] to compute $\langle S_{2}\rangle$ and $\langle V_{2}\rangle$ exactly for $N\geq 1$ planar Brownian motions. Since then, it has been used for a number of continuous-time planar processes: random acceleration process [17], branching Brownian motion with applications to animal epidemic outbreak [18], anomalous diffusion processes [19] and also to a single Brownian motion confined to a half-space [20].

The objective of this paper is to go beyond the continuous-time limit and obtain results for the convex hull of a discrete-time planar random walk of $n$ steps (with $n$ large but finite) with arbitrary jump distribution, including for instance Lévy flights. Indeed, in any realistic experiment or simulation, the points of the trajectory are always discrete. For example, recently proposed local convex hull estimators [21] are based on a relatively small number of points, where we can not apply the Brownian limiting results reviewed above. The first rigorous result for a two-dimensional discrete random walk, modeled as a sum of independent random variables in the complex plane, was derived for the mean perimeter of the convex hull by Spitzer and Widom [22],

[TABLE]

(here $x_{k}+iy_{k}$ is the complex-valued position of the walker after $k$ steps). Although the formula (3) looks deceptively simple, an explicit computation of the mean $\langle L_{n}\rangle$ is difficult using Eq. (3), in particular its behavior for large but finite $n$. For the case of zero mean and finite variance jump distributions, the leading $n^{\frac{1}{2}}$ term in $\langle L_{n}\rangle$ was identified in [22], but the relevant subleading terms were not known, to the best of our knowledge. For other results on the statistics of $L_{n}$, see Ref. [23]. The Spitzer-Widom formula (3) was extended to generic $d$-dimensional random walks in Ref. [24], in which exact combinatorial expressions for the expected surface area and the expected volume of the convex hull were derived. However, these expressions are not suitable for the asymptotic analysis at large $n$. Several other geometrical properties of the convex hull of random walks are known: for example, the exact formula for the mean number of facets of the convex polytope of a $d$-dimensional random walk, for $d=2$ [25] and $d>2$ [26, 27]. But in this paper, we will restrict ourselves only to the mean perimeter $\langle L_{n}\rangle$ and the mean area $\langle A_{n}\rangle$ of a planar random walk of $n$ steps and our main goal is to derive explicitly not only the leading term in $\langle L_{n}\rangle$ and $\langle A_{n}\rangle$ for large $n$, but also the subleading terms.

Our strategy is to adapt the mapping between the convex hull of a $2$-d process and the extreme statistics of the $1$-d component process, mentioned above, to the case of a single discrete-time planar random walk with generic jump distributions. Using this strategy, we are able to compute explicitly the leading and subleading terms of the mean perimeter of the convex hull of a planar random walk of $n$ steps with arbitrary symmetric continuous jump distributions for large but finite $n$. The mean area is also computed but only for the particular case of isotropic Gaussian jumps. The rest of the paper is organized as follows. Section 2 outlines the class of considered planar random walks and the main results. In Sec. 3, we explain the derivation steps. In Sec. 4, several explicit examples are presented and used to illustrate the accuracy of the derived asymptotic relations by comparison with Monte Carlo simulations. In Sec. 5, we discuss the main results, their applications, and conclusions. A and B contain some technical details of the derivation and exactly solvable examples, respectively.

2 The model and the main results

We consider a discrete-time random walker in the plane whose jumps are random, independent, and identically distributed. Starting from the origin, the walker produces a sequence of $(n+1)$ points $\{(x_{0},y_{0}),(x_{1},y_{1}),\ldots,(x_{n},y_{n})\}\subset\mathbb{R}^{2}$ after $n$ jumps such that

[TABLE]

where the jumps $\vec{\eta}_{k}=(\eta_{k}^{x},\eta_{k}^{y})$ at the $k$-th step are independent from step to step, and at each step they are drawn from a prescribed joint probability density function (PDF) $p(x,y)$, i.e.,

[TABLE]

We emphasize that the starting point $(x_{0},y_{0})$ is not random and for convenience, we choose $(x_{0}=y_{0}=0)$ to be the origin. The convex hull constructed over these $(n+1)$ points is the minimal convex polygon that encloses all these points (see Fig. 1 for an illustration). We are interested in the perimeter $L_{n}$ and the area $A_{n}$ of the convex hull which are random variables given that the points are generated as successive positions of a planar random walk. We aim at computing exactly the leading and subleading terms of the mean perimeter, $\langle L_{n}\rangle$, and the mean area, $\langle A_{n}\rangle$, of the convex hull for large $n$.

As explained in Sec. 3, our computation relies on two key results: (i) Cauchy’s formula for the perimeter and the area of a closed convex curve, that allows one to reduce the original planar problem to the analysis of one-dimensional projections, and (ii) the Pollaczek-Spitzer formula describing the distribution of the maximum of partial sums of independent symmetric continuously distributed random variables [28, 29]. To use the Pollaczek-Spitzer formula, we need thus to assume that the joint probability density $p(x,y)$ is continuous and centrally symmetric:

[TABLE]

In particular, our results will not be applicable to a classical random walk on the square lattice because its distribution is not continuous. In the following, we outline the main results that will be derived in Sec. 3.

The mean perimeter $\langle L_{n}\rangle$ is computed for a very general class of symmetric continuous jump distributions. Writing the Fourier transform of $p(x,y)$ as

[TABLE]

one can characterize the behavior of the mean perimeter according to the asymptotic properties of $\hat{\rho}_{\theta}(k)$ as $k\to 0$. We assume a general expansion

[TABLE]

with the scaling exponent $0<\mu\leq 2$, and a scale $a_{\theta}>0$. When $0<\mu\leq 1$, the mean perimeter of the convex hull is infinite. We therefore focus on the case $1<\mu\leq 2$.

First, we derive an exact formula for the generating function of $\langle L_{n}\rangle$ which is valid for any $1<\mu\leq 2$. Extracting the asymptotic large $n$ behavior of $\langle L_{n}\rangle$ from this general formula is, however, nontrivial. We distinguish the following cases.

When the jump variance is finite ($\mu=2$), the mean perimeter is shown to behave as

[TABLE]

with

[TABLE]

where $\sigma_{\theta}$ and $\gamma_{\theta}$ are given by

[TABLE]

If in addition the fourth-order moment of the jump distribution is finite, one gains the second subleading term,

[TABLE]

with

[TABLE]

and

[TABLE]

Higher-order corrections can also be derived under further moments assumptions. Note that the integral expression for the coefficient $C_{0}$ in front of the leading term $n^{1/2}$ first appeared in [22]. In Sec. 4, we will show that the asymptotic formula (13) is very accurate even for small $n$. 2. 2.

When the jump variance is infinite (i.e., $1<\mu<2$), one needs to consider the subleading term in the small $k$ asymptotics of $\hat{\rho}_{\theta}(k)$:

[TABLE]

with the subleading exponent $\nu>\mu$ and a coefficient $b_{\theta}$. Depending on the subleading exponent $\nu$, we distinguish two cases:

(1) if $\mu<\nu<\mu+1$, one has

[TABLE]

with

[TABLE]

Note that the coefficient $C_{0}$ also appears in the mean perimeter of the convex hull of continuous-time symmetric stable processes [13].

(2) if $\nu>\mu+1$, one has

[TABLE]

with $C_{0}$ from Eq. (18) and

[TABLE]

where

[TABLE]

For instance, for a Lévy symmetric alpha-stable distribution with $\hat{\rho}(k)=\exp(-|ak|^{\mu})$, one gets [33]

[TABLE]

where $\zeta(z)$ is the Riemann zeta function.

The obtained results are indeed very general.

In turn, our method of computation of the mean area requires two additional strong assumptions: (a) the independence of the jumps along $x$ and $y$ coordinates, i.e., $p(x,y)=p(x)p(y)$ and (b) the isotropy of the jump PDF, i.e., $p(x,y)$ should depend only on the distance $r=\sqrt{x^{2}+y^{2}}$ but not on the direction of the jump. According to Porter-Rosenzweig theorem [30], only the Gaussian jump distribution with identical variance $\sigma^{2}$ along $x$ and $y$ directions, i.e., $p(x,y)=\frac{1}{2\pi\sigma^{2}}\,\exp[-(x^{2}+y^{2})/{2\sigma^{2}}]$, satisfies both properties (a) and (b). Our result for the mean area is thus only valid for this Gaussian distribution:

[TABLE]

with ${\mathcal{K}}=3$ and

[TABLE]

Recently, an exact formula the mean area of the convex hull of a Gaussian random walk was derived [8]. In the isotropic case, the formula reads

[TABLE]

While the result in Eq. (26) is very useful for finite $n$, deriving the large $n$ asymptotics of this double sum (including up to two subleading terms as in Eq. (24)) seems somewhat complicated. Our method, in contrast, gives a more direct access to the asymptotics. Moreover, one can check numerically that our asymptotic formula (24) agrees accurately with the exact expression (26) even for moderately large $n$.

The leading term of Eq. (24) was shown to be valid for a generic random walk with increments of a finite variance (see Proposition 3.3 in [9]). Moreover, our numerical simulations (see Sec. 4.2 and Fig. 5) suggest that the obtained formula (24) (including the subleading terms) may be applicable for some other isotropic processes. In other words, the technical assumption about the independence of the jumps along $x$ and $y$ might be relaxed in future. This statement, which is uniquely based on numerical simulations for some jump distributions, is conjectural. In turn, the isotropy assumption is important, as illustrated by numerical simulations.

The large $n$ asymptotic relations (13, 17, 20, 24) are the main results of the paper. Setting $t=n\tau$ and $D=2\sigma^{2}/(4\tau)$ with a time step $\tau$, one recovers from Eqs. (13, 24) the same leading terms as in Eq. (1, 2) for Brownian motion (note that we write $2\sigma^{2}$ in $D$ because $\sigma^{2}$ is the variance of jumps along one direction). It is thus not surprising that the leading term in Eq. (13) is universal because its derivation is valid for any planar random walk with a symmetric and continuous jump distribution and having a finite variance $\sigma^{2}$.

Thinking of Brownian motion as a limit of random walks, the subleading terms in Eq. (13) can be understood as “finite size” corrections. The first subleading term is valid under the same assumptions as the leading term, although the coefficient $\gamma_{\theta}$ depends on the jump distribution (see examples in Sec. 4). In turn, the second subleading term depends on the kurtosis ${\mathcal{K}}_{\theta}$ and thus requires an additional assumption that ${\mathcal{K}}_{\theta}$ is finite.

3 Main steps leading to the derivation of results

3.1 Reduction to a one-dimensional problem

We start with Cauchy’s formula for the perimeter $L$ and the area $A$ of an arbitrary convex domain ${\cal C}$ with a reasonably smooth boundary $\gamma_{\cal C}$ [15, 16]. Let the boundary $\gamma_{\cal C}$ be parameterized as $(X(s),Y(s))$ with a curvilinear coordinate $s$ ranging from [math] to $1$. Setting the origin of coordinates inside the domain, one defines the support function $M(\theta)$ as the distance from the origin to the closest straight line that does not cross the domain and is perpendicular to the vector from the origin in direction $\theta$ (Fig. 2). In other words,

[TABLE]

Cauchy showed that [15]

[TABLE]

For a simple derivation of this formula see Ref. [1]. A straightforward calculation of $M(\theta)$ for a convex hull over a set of points may seem to be hopeless, as one would need first to construct the convex hull by identifying and ordering its vertices among the given set of points and then to compute $M(\theta)$. The key idea is that $M(\theta)$ can be found directly from the vertices of the trajectory as [22, 6]

[TABLE]

Moreover, given that the maximum for a fixed $\theta$ is realized by a certain vertex (with index $k^{*}$ which discretely changes with $\theta$), one also obtains the derivative:

[TABLE]

When the points $(x_{k},y_{k})$ are random, the perimeter and the area of the convex hull are random variables. We focus on the mean values $\langle L_{n}\rangle$ and $\langle A_{n}\rangle$:

[TABLE]

i.e., the computation is reduced to the first two moments of $M(\theta)$ and to the mean $\langle[M^{\prime}(\theta)]^{2}\rangle$. The important observation is that, for a fixed direction $\theta$, one needs to characterize the maximum of the projection of points $(x_{k},y_{k})$ onto that direction

[TABLE]

The projection of a random walk is also a random walk. In fact, we can write according to Eq. (4)

[TABLE]

with

[TABLE]

The probability density of $\xi_{k}^{\theta}$, $\rho_{\theta}(z)$, is fully determined by that of the jump $(\eta_{k}^{x},\eta_{k}^{y})$. In particular, its Fourier transform $\hat{\rho}_{\theta}(k)$ is related to $p(x,y)$ by Eq. (7). The symmetry (6) implies that $\hat{\rho}_{\theta}(-k)=\hat{\rho}_{\theta}(k)$ and thus the density $\rho_{\theta}(z)$ is symmetric.

Having discussed the general jump distributions, let us mention two particular cases that will be important later.

(a) in the case of independent jumps along $x$ and $y$ coordinates, one has $p(x,y)=p_{x}(x)p_{y}(y)$, and thus

[TABLE]

where $\hat{p}_{x}$ and $\hat{p}_{y}$ are the Fourier transforms of $p_{x}(x)$ and $p_{y}(y)$, respectively.

(b) For isotropic jumps, $p(x,y)$ depends only on the radial coordinate, $p(x,y)dxdy=p_{r}(r)dr\,d\phi/(2\pi)$, where $p_{r}(r)$ is the radial density (that includes the factor $r$ from the Jacobian). From Eq. (7), one gets

[TABLE]

in which the integration over the angular coordinate $\phi$ eliminated the dependence on $\theta$ and resulted in the Bessel function of the first kind, $J_{0}(|k|r)$.

3.2 Formal solution of the one-dimensional problem

The formal exact solution of the one-dimensional problem can be obtained via the Pollaczek-Spitzer formula [28, 29]. This formula characterizes the maximum of partial sums of independent identically distributed random variables $\xi_{k}$ with a symmetric and continuous density $\rho(z)$:

[TABLE]

(in this subsection, we temporarily drop the subscript and superscript $\theta$ from all the variables; it will be restored at the end). Considering $\xi_{k}$ as jumps of a random walker, $z_{k}=z_{k-1}+\xi_{k}$ (with $z_{0}=0$), one can also write

[TABLE]

Pollaczek and later Spitzer showed that the cumulative distribution $Q_{n}(z)=\P\{M_{n}\leq z\}$ of $M_{n}$ satisfies the following identity for $0\leq s\leq 1$ and $\lambda\geq 0$

[TABLE]

with

[TABLE]

where $Q^{\prime}_{n}(z)=dQ_{n}(z)/dz$ is the probability density of the maximum [28, 29]. In principle, all the moments of $M_{n}$ can be obtained from the formula in Eq. (41). However, in practice, deriving explicitly the moments of $M_{n}$ by inverting this formula is highly nontrivial [31]. For example, the expected maximum of a discrete-time random walk $\langle M_{n}\rangle$ appears in a number of different problems, from packing algorithms in computer science [32], all the way to the survival probability of a single or multiple walkers in presence of a trap [33, 34, 35, 36, 37, 38]. It has also appeared in the context of the order, gap and record statistics of random walks [39, 40, 41, 42, 43] and $\langle M_{n}\rangle$ has been analyzed for large $n$ (for the leading and the next subleading term) in detail using the Pollaczek-Spitzer formula in Eq. (41). Here, in addition to calculating the first three terms in the asymptotic expansion of $\langle M_{n}\rangle$ for $n\gg 1$, we also calculate the large $n$ behavior of the second moment $\langle M_{n}^{2}\rangle$, that we need for the computation of the mean area of the convex hull.

In fact, the Pollaczek-Spitzer formula also determines the generating functions for all moments of $M_{n}$:

[TABLE]

By considering a general asymptotic expansion

[TABLE]

with an exponent $0<\mu\leq 2$ and a scale $a>0$, we derive in A.1 the exact expressions

[TABLE]

which is valid for any $1<\mu\leq 2$ (note that $\langle M_{n}\rangle=\infty$ for $0<\mu\leq 1$), and

[TABLE]

which is valid for $\mu=2$ (note that $\langle M_{n}^{2}\rangle=\infty$ for $0<\mu<2$). The exact relations (45, 46) are new results which allow one to study the first two moments of the maximum $M_{n}$. In the next subsection, we will analyze the expansion of Eqs. (45, 46) as $s\to 1$ in order to determine the asymptotic behavior of the moments $\langle M_{n}\rangle$ and $\langle M_{n}^{2}\rangle$ as $n\to\infty$. We consider separately jumps with a finite variance, and Lévy flights.

3.3 Mean perimeter and mean area of the convex hull

3.3.1 Mean perimeter for jumps with a finite variance.

Given a generic continuous jump distribution $p(x,y)$ satisfying the property in Eq. (6), we determine $\hat{\rho}_{\theta}(k)$ using Eq. (7). Furthermore, by examining the small $k$ behavior of $\hat{\rho}_{\theta}(k)$, we determine the $\theta$-dependent variance $\sigma_{\theta}^{2}$ and the $\theta$-dependent kurtosis ${\mathcal{K}}_{\theta}$, using respectively Eqs. (11) and (15). In addition, knowing $\hat{\rho}_{\theta}(k)$ from Eq. (7), we also determine $\gamma_{\theta}$ in Eq. (12). Equipped with these three quantities $\sigma_{\theta}$, ${\mathcal{K}}_{\theta}$ and $\gamma_{\theta}$, we show in A that the leading large $n$ terms of the first two moments of $M_{n}$ are given by

[TABLE]

For the mean perimeter of the convex hull, we will only need the first moment in Eq. (47). Indeed, using Eq. (32), the integration of the expansion (47) over $\theta$ from [math] to $2\pi$ yields the announced result (13) for the mean perimeter of the convex hull. The result for the second moment in Eq. (48) will be needed later to determine the mean area $\langle A_{n}\rangle$.

3.3.2 Mean perimeter for Lévy flights.

When the variance of jumps is infinite, one gets the Taylor expansion Eq. (8), with the scaling exponent $0<\mu<2$. When $0<\mu\leq 1$, the mean perimeter of the convex hull is infinite. Throughout this section, we focus on the case $1<\mu<2$, in which the first moment of jumps is finite (and zero due to the assumption of a symmetric distribution), whereas the variance is infinite. In this case, the leading behavior of the mean maximum of partial sums is universal [33]

[TABLE]

However, the subleading term depends on finer details of the jump distribution. In order to determine the subleading term, we consider the expansion (16) with the subleading term $b_{\theta}|k|^{\nu}$ such that $\nu>\mu$. We distinguish two cases: $\mu<\nu<\mu+1$ and $\nu>\mu+1$. In A.3, we derive the following asymptotics results:

[TABLE]

for $\mu<\nu<\mu+1$, and

[TABLE]

for $\nu>\mu+1$, with $\gamma_{\theta}$ given by Eq. (22). The asymptotic relation (51) was first derived in [33] for the particular case $\nu=2\mu$. One can see that for $\mu<\nu<\mu+1$, the subleading term of $\langle M_{n}\rangle$ grows with $n$, whereas for when $\nu>\mu+1$, the subleading term is constant. Higher-order corrections can be derived as well.

Finally, using the Cauchy formula (32), the integration of Eqs. (50, 51) over $\theta$ from [math] to $2\pi$ yields Eqs. (17, 20) for the mean perimeter of the convex hull, announced in Section 2.

3.3.3 Mean area for isotropic Gaussian jumps.

According to Eq. (33), the expansion (48) determines the first contribution to the mean area. This contribution was calculated for an arbitrary symmetric continuous jump distribution with a finite variance. The second contribution to the mean area comes from $\langle[M^{\prime}(\theta)]^{2}\rangle$ that has to be computed separately. We recall that $M^{\prime}(\theta)$ is given by Eq. (31). Our computation of this contribution relies on two additional simplifying assumptions: (a) the jumps along $x$ and $y$ coordinates are independent and (b) the jump process is isotropic, i.e., the distribution of jumps does not depend on their direction. In this case, using the isotropy condition (b), we get

[TABLE]

The disentanglement of $\langle[M(\theta)]^{2}\rangle$ and $\langle[M^{\prime}(\theta)]^{2}\rangle$ allows one to compute the latter one by using the following argument. We recall that $k^{*}$ is the index of the maximal position among $x_{k}$, i.e., its statistics is fully determined by the jumps along $x$ coordinate. Once this statistics is known, the mean $\langle[M^{\prime}(0)]^{2}\rangle$ can be found by taking the conditional expectation of $y^{2}_{k^{*}}$ at any fixed value $k^{*}$ and then the expectation with respect to the distribution of $k^{*}$. Now, once we condition $k^{*}$, i.e., the time step at which the $x_{k}$’s achieve their maximum, the $y_{k}$ process will, in general, be affected by this conditioning. However, if $x_{k}$ and $y_{k}$ are independent (which happens when the jump process satisfies property (a) above), we get

[TABLE]

where $\langle[\eta_{k}^{y}]^{2}\rangle=\sigma^{2}$. It remains to find $\langle k^{*}\rangle$. For symmetric and continuous jump distributions, it follows from symmetry that $\langle k^{*}\rangle=n/2$ independent of the details of the jump PDF. This can be deduced formally also, by noting that the time step $k^{*}$, at which the maximum of $x_{k}$ is achieved, has a universal distribution, independent of the jump distribution (given that the latter is continuous and symmetric) [31]:

[TABLE]

This is the direct consequence of the Sparre Andersen theorem. From this distribution, one easily computes the mean value as

[TABLE]

and thus

[TABLE]

This yields Eq. (24) for the mean area of the convex hull for an isotropic jump process with independent jumps along $x$ and $y$ coordinates. As discussed in Sec. 2, the only process satisfying two requirements (a) and (b) is the isotropic Gaussian process. Our formula (24) is thus provided only in this case. For a more general case, one would need to compute the joint distribution of the maximum position $k^{*}$ and both values $x_{k^{*}}$ and $y_{k^{*}}$ which is more complicated and remains an open problem.

4 Examples and simulations

In this section, we illustrate the above general results on several examples of symmetric planar random walks. We derive the explicit values of the relevant parameters that determine the mean perimeter and the mean area. We also investigate the effect of anisotropy of the jump distribution on the convex hull properties. Finally, we compare our theoretical predictions to the results of Monte Carlo simulations. For this purpose, we generate $10^{5}$ planar random walks, compute the convex hull for each generated trajectory with $n+1$ points by using the Matlab function ’convhull’, and determine its perimeter and area. These simulations yield the representative statistics of perimeters and areas from which the mean values are computed.

4.1 Gaussian jumps

We first consider the basic example of Gaussian jumps which are independent along $x$ and $y$ coordinates and characterized by variances $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$. Substituting the jump probability density,

[TABLE]

into Eq. (7) yields $\hat{\rho}_{\theta}(k)=e^{-k^{2}\sigma^{2}_{\theta}/2}$, with $\sigma_{\theta}^{2}=\sigma_{x}^{2}\cos^{2}\theta+\sigma_{y}^{2}\sin^{2}\theta$. One can see that the anisotropy only affects the variance $\sigma_{\theta}^{2}$, whereas the two other relevant parameters, $\gamma$ and ${\mathcal{K}}$, which are rescaled by variance, do not depend on $\theta$. One finds ${\mathcal{K}}=3$, whereas the integral in Eq. (12) was computed exactly in [33] and provided in Eq. (25). Assuming (without loss of generality) that $\sigma_{x}\geq\sigma_{y}$, we set

[TABLE]

where $E(k)$ is the complete elliptic integral of the second kind (for the isotropic case, $\sigma_{x}=\sigma_{y}=\sigma$). With this notation, we get the expansion coefficients

[TABLE]

Figure 3 shows the rescaled mean perimeter, $\langle L_{n}\rangle/n^{1/2}$, and the rescaled mean area, $\langle A_{n}\rangle/n$, for isotropic planar random walk with independent Gaussian jumps ($\sigma_{x}=\sigma_{y}=1$). The results of Monte Carlo simulations are in perfect agreement with our theoretical predictions (13, 24). One can see that the subleading terms play an important role. In fact, if one kept only the leading term and omitted the subleading terms, one would get the horizontal dotted line. This corresponds to the case of Brownian motion, in which only the leading term is present, see Eqs. (1, 2). The subleading terms account for the discrete-time character of random walks which is particularly important for moderate values of $n$. In order to highlight the role of the third term in the asymptotic expansions, we also draw by dashed line the asymptotics without this term. As expected, the third term improves the quality of the theoretical prediction at small $n$. Note also that the asymptotic relation is slightly less accurate for the mean area.

In Fig. 4, we consider the convex hull for anisotropic random walk with independent Gaussian jumps, with $\sigma_{x}=5$ and $\sigma_{y}=1$. In this case, one can use the asymptotic formula (13), in which the expansion coefficients $C_{j}$ are given by Eq. (60), with an effective variance $\sigma^{2}$ computed in Eq. (59). In this example, \sigma=5\frac{2}{\pi}E\bigl{(}\sqrt{1-1/25}\bigr{)}=3.3439\ldots. For the mean perimeter, one observes a perfect agreement between the theoretical predictions and Monte Carlo simulations. In turn, our asymptotic formula (24) for the mean area is not applicable for anisotropic case, as also confirmed by simulations (not shown).

4.2 Exponentially distributed radial jumps

The next common model has exponentially distributed radial jumps with uniform angular distribution. This is a particular realization of a “run-and-tumble” model of bacterial motion [44, 45, 46]. Substituting the radial density $p_{r}(r)=\sigma^{-1}\,e^{-r/\sigma}$ into Eq. (38) yields \hat{\rho}(k)=\bigl{(}1+(k\sigma)^{2}\bigr{)}^{-1/2}. One gets ${\mathcal{K}}=9$ and

[TABLE]

from which the expansion coefficients are

[TABLE]

Figure 5 illustrates the obtained results. As for isotropic Gaussian jumps, there is a perfect agreement between theory and simulations for the mean perimeter. We also present the results for the mean area. We recall that the theoretical formula (24) was derived under the assumption of independent jumps along $x$ and $y$ coordinates, which evidently fails for the exponential radial jumps. In spite of this failure, the theoretical formula (24) is in perfect agreement with Monte Carlo simulations, except for very small $n$. This empirical observation suggests a possibility to relax this technical assumption in future, at least for large $n$.

4.3 Independent exponentially distributed jumps

We also consider the case, when the jumps along $x$ and $y$ coordinates are independent and exponentially distributed, with two densities $p_{x}(x)=\frac{1}{2}\sigma_{x}^{-1}e^{-|x|/\sigma_{x}}$ and $p_{y}(y)=\frac{1}{2}\sigma_{y}^{-1}e^{-|y|/\sigma_{y}}$. The Fourier transforms are $\hat{p}_{x}(k)=(1+(k\sigma_{x})^{2})^{-1}$ and $\hat{p}_{y}(k)=(1+(k\sigma_{y})^{2})^{-1}$ so that

[TABLE]

We get $\sigma^{2}_{\theta}=2(\sigma_{x}^{2}\cos^{2}\theta+\sigma_{y}^{2}\sin^{2}\theta)$ and

[TABLE]

Using the identity (111), we compute explicitly

[TABLE]

In the particular case $\sigma_{x}=\sigma_{y}=\sigma$, one gets $\sigma^{2}_{\theta}=2\sigma^{2}$, ${\mathcal{K}}_{\theta}=6(1-\sin^{2}\theta\cos^{2}\theta)$, and

[TABLE]

For this case, we obtain from Eqs. (10)

[TABLE]

Figure 6 illustrates the excellent agreement between theory and simulations.

4.4 Radial Lévy jumps

Now, we investigate the example of radial Lévy flights with infinite variance (but finite mean) and uniform angle distribution. Among various heavy-tailed jump distributions (e.g., Pareto distributions), we choose for our illustrative purposes the distribution

[TABLE]

with a scale $R>0$ and the scaling exponent $\mu=2\alpha$, with $\frac{1}{2}<\alpha<1$. For this distribution, Eq. (38) yields a simple closed formula [47]

[TABLE]

where $K_{\alpha}(z)$ is the modified Bessel function of the second kind. The asymptotic behavior of $K_{\alpha}(z)$ near $z$ implies as $k\to 0$

[TABLE]

Comparing this expansion to Eq. (16), one can identify

[TABLE]

In Fig. 7, the mean perimeter computed by Monte Carlo simulations for $\mu=1.5$ and $R=1$ is compared to the theoretical prediction (17). One can see that the agreement is good but worse than for the earlier examples with a finite variance. One obvious reason is that here we have determined only two terms, whereas Eq. (13) contains three terms.

4.5 Independent Lévy $\alpha$-stable symmetric jumps

Finally, we investigate Lévy $\alpha$-stable symmetric jumps, with independent displacements along $x$ and $y$ coordinates given by $\hat{p}_{x}(k)=\hat{p}_{y}(k)=\exp(-|ak|^{\mu})$, with $1<\mu<2$. Using Eq. (37), one gets thus $\hat{\rho}_{\theta}(k)=\exp(-|a_{\theta}k|^{\mu})$, with

[TABLE]

so that $\nu=2\mu$, and $b_{\theta}=a_{\theta}^{2}/2$. The mean perimeter is determined by Eq. (20), with

[TABLE]

and $\gamma$ given by Eq. (23) which is independent of $\theta$.

For $\mu=3/2$, we obtain numerically $a^{-1}C_{0}=8.6275\ldots$ and $a^{-1}C_{1}=-5.2151\ldots$. Figure 8 shows the good agreement between the theoretical prediction (20) and Monte Carlo simulations.

5 Discussion and conclusion

To summarize, we have presented exact asymptotic results for the mean perimeter of the convex hull of an $n$-step discrete-time random walk in a plane, with a generic continuous jump distribution satisfying the central symmetry assumption in Eq. (6). Explicit results, along with simulations confirming them have been presented for several examples of such jump distributions. For the mean area of the convex hull, we have derived exact results for isotropic Gaussian jump distributions. For jumps with a finite variance, our results provide precise estimates of the deviations from the Brownian limit and explain the discrepancies between the asymptotic Brownian limit results and observed simulations for finite but large $n$.

The obtained results are particularly valuable for applications dealing with discrete-time random processes, e.g., home range estimation in ecology. Given that the tracks of animal displacements are typically recorded at discrete time steps (e.g., daily observations) and relatively short, the subleading terms play an important role. The asymptotic formulas can also be used for calibrating new estimators, based on the local convex hull, that were proposed for the analysis of intermittent processes in microbiology [21]. Finally, the knowledge of the mean perimeter of the convex hull can be used to estimate the scaling exponent and the scale of symmetric Lévy flights, for which the conventional mean and variance estimators are useless.

There are many interesting open problems that may possibly be addressed using the methods presented here. For example, the numerical evidence suggests a possible extension of the derived asymptotic formula for the mean area to other isotropic processes, beyond the Gaussian case. Also, it would be interesting to extend our results to the case of the convex hull of planar discrete-time random bridges (where the walker is constrained to come back to the starting point after $n$ discrete jumps). For such discrete-time bridges, there are recent exact results on the statistics of the first two maximum and the gap between them [48] which may be useful for the convex hull problem. One can also consider the problem with many independent discrete-time walkers. Finally, it would be interesting to study the statistics of the perimeter and the area for random walks with jump distributions that violate the reflection property in Eq. (6), for example, for walks in presence of a drift or a potential.

Acknowledgments

DG acknowledges the support under Grant No. ANR-13-JSV5-0006-01 of the French National Research Agency.

Appendix A Asymptotic analysis

In this Appendix, we derive the main results (47, 48). The derivation is based on the asymptotic analysis of the Pollaczek-Spitzer identity (41) and extends the earlier results from [33]. Since the function $\phi(s,\lambda)$ from Eq. (42) is related to the Laplace transform of the probability density $Q^{\prime}_{n}(z)$, it determines the generating functions $h_{m}(s)$ of the moments $\langle M_{n}^{m}\rangle$ via Eq. (43). In turn, the asymptotic behavior of $h_{m}(s)$ as $s\to 1$ determines the asymptotic behavior of $\langle M_{n}^{m}\rangle$ as $n\to\infty$.

A.1 Derivation of generating functions

We recall that $M_{n}$ denotes the maximum of partial sums of independent identically distributed random variables $\{\xi_{k}\}$:

[TABLE]

We assume that the jump distribution is symmetric and continuous, whereas its characteristic function $\hat{\rho}(k)$ admits the expansion

[TABLE]

with an exponent $0<\mu\leq 2$ and a scale $a>0$.

First, we derive the generating function $h_{1}(s)$ from the Pollaczek-Spitzer formula (41). For this purpose, we need to determine the expansion of this formula in powers of $\lambda$ for small $\lambda$. Let us first write

[TABLE]

It is easy to obtain the $\lambda\to 0$ limit of $I(s,\lambda)$. We rescale $k=\lambda\,u$ in Eq. (76) and take the $\lambda\to 0$ limit. This gives, using ${\hat{\rho}}(k=0)=1$, a very simple expression

[TABLE]

Next we re-write

[TABLE]

Expanding the right-hand side of Eq. (78) up to the first order $\lambda$ for small $\lambda$ (with $s$ fixed) gives

[TABLE]

Note that the integral in the second term is convergent for any jump pdf with ${\hat{\rho}}(k)$ satisfying Eq. (75) provided $1<\mu\leq 2$. Taking the derivative with respect to $\lambda$, evaluated at $\lambda=0$, yields

[TABLE]

Note that this is an exact result for any jump PDF ${\hat{\rho}}(k)$ satisfying Eq. (75) with $1<\mu\leq 2$ and it also holds for arbitrary $s$ such that the generating function on the left-hand side of Eq. (80) is convergent. In turn, when $0<\mu\leq 1$, the mean value $\langle M_{n}\rangle$ is infinite, and $h_{1}(s)$ is undefined.

The second moment $\langle M_{n}^{2}\rangle$ is finite only for $\mu=2$. Expanding Eq. (78) up to the order $\lambda^{2}$ for small $\lambda$ (with $s$ fixed) and taking the second derivative with respect to $\lambda$ at $\lambda=0$, one gets for any $0\leq s<1$

[TABLE]

A.2 Jumps with a finite variance

Here we investigate the asymptotic behavior of $h_{1}(s)$ as $s\to 1$ for the jump distribution with a finite variance $\sigma^{2}$ so that $\hat{\rho}(k)\simeq 1-k^{2}\sigma^{2}/2+o(k^{2})$. First, we represent $h_{1}(s)$ from Eq. (45) as

[TABLE]

The first integral can be computed explicitly, so that

[TABLE]

where $I_{\varepsilon}$ denotes the second integral in Eq. (82) that we rewrite as

[TABLE]

with $\varepsilon=\sqrt{1-s}$. Setting $\varepsilon\to 0$, one gets

[TABLE]

with the constant $\gamma$ from Eq. (12), with $\sigma=\sigma_{\theta}$.

Next, we consider

[TABLE]

Changing again the integration variable, one has

[TABLE]

To get the next-order term, we assume the existence of the fourth-order moment of the jumps so that

[TABLE]

where $c={\mathcal{K}}/24$ is related to the kurtosis ${\mathcal{K}}$ given by Eq. (15). Substituting this expansion into Eq. (87), one gets

[TABLE]

Combining these results, we conclude that

[TABLE]

Substituting this expansion into Eq. (46), one gets

[TABLE]

In order to invert the relations (90, 91), we use the following identities:

[TABLE]

with

[TABLE]

Inverting term by term, we obtain the main relations (47, 48). Higher-order corrections can be obtained in the same way.

A.3 Lévy flights

For Lévy flights with $\hat{\rho}(k)\simeq 1-|ak|^{\mu}+o(|k|^{\mu})$, one can derive the asymptotic behavior of $h_{1}(s)$ as $s\to 1$ for the case $1<\mu<2$. In this case, we replace the representation (82) by

[TABLE]

As previously, the first integral can be evaluated explicitly so that

[TABLE]

where

[TABLE]

with $\varepsilon=\sqrt{1-s}$. The first term in Eq. (101) provides the leading contribution. Using the discrete form of the Tauberian theorem,

[TABLE]

one gets the leading term to $\langle M_{n}\rangle$ to be $(a\mu\Gamma(1-1/\mu)/\pi)\,n^{1/\mu}$.

In order to compute the subleading term, one needs to analyze the integral $I_{\varepsilon}$ whose behavior depends on the next-order term in the small $k$ expansion of $\hat{\rho}(k)$ in Eq. (16). Formally setting $\varepsilon=0$ in Eq. (102), one would get

[TABLE]

Substituting the expansion (16) into this integral, one gets the factor $\ln(1-bk^{\nu-\mu}/a^{\mu})$. As a consequence, the integral converges when $\nu-\mu>1$ and diverges otherwise. This condition naturally distinguishes two asymptotic regimes: $\nu>\mu+1$ and $\nu\leq\mu+1$. In the former case, the integral converges, and the two principal terms in $h_{1}(s)$ are

[TABLE]

where $I_{0}$ is identical to $\gamma$ from Eq. (22). The inversion of this relation yields Eq. (51).

In the case $\nu\leq\mu+1$, the subtraction of the term with $1-|ak|^{\mu}$ in Eq. (100) was not sufficient to get the convergent correction term. We remedy this problem by using another representation:

[TABLE]

Changing the variable $u=ak(s/(1-s))^{1/\mu}$, we rewrite the first term as

[TABLE]

with $C=ba^{-\nu}s^{-\nu/\mu}(1-s)^{\nu/\mu}$. When $s\to 1$, $C$ is small, and one can expand the logarithm to get

[TABLE]

Now both integrals are convergent. The first integral yields the same leading term as earlier, whereas the second term provides the subleading correction. Finally, the second term in Eq. (106) results in higher-order corrections that we ignore. Keeping only the leading term and the first subleading term, we get as $s\to 1$

[TABLE]

The inversion of this relation yields the large $n$ asymptotic formula (50) for $\langle M_{n}\rangle$.

Appendix B Two exactly solvable examples

In this Appendix, we briefly discuss two distributions, for which the function $\phi(s,\lambda)$ can be found exactly: the symmetric exponential distribution $\rho(z)=\frac{1}{2}e^{-|z|}$ (see [33]) and its modification $\rho(z)=\frac{1}{2}|z|e^{-|z|}$. These examples served us as benchmarks for checking general results.

B.1 Symmetric exponential distribution

When $\rho(z)=\frac{1}{2}e^{-|z|}$, one gets $\hat{\rho}(k)=(1+k^{2})^{-1}$ and thus Eq. (42) becomes

[TABLE]

where we used the identity

[TABLE]

to compute the integral:

[TABLE]

According to Eqs. (11, 15, 12), we also obtain

[TABLE]

The first derivative of $\phi(s,\lambda)$ from Eq. (110) yields

[TABLE]

from which, using Eqs. (92, 93), one retrieves the exact form of the first moment, which is valid for any $n$ and was first derived in [33]:

[TABLE]

Using the Stirling formula, one can get the large $n$ expansion to any order:

[TABLE]

The first three terms agree with the general expansion (47).

Similarly, we get from Eq. (110)

[TABLE]

from which the exact formula follows using Eqs. (92, 93)

[TABLE]

In the large $n$ limit, we deduce

[TABLE]

The first three terms agree with the general expansion (48).

Note that one can also obtain the generating function for the cumulative distribution $Q_{n}(z)$. In fact, one has for the Laplace transform of $Q_{n}(z)$

[TABLE]

from which the inverse Laplace transform (with respect to $\lambda$) yields

[TABLE]

From this relation, one can easily get any generating function $h_{m}(s)$.

B.2 Modified exponential distribution

For the case $\rho(z)=\frac{1}{2}|z|e^{-|z|}$, one gets

[TABLE]

Using again Eq. (111), one computes the integral in Eq.(42) as

[TABLE]

where

[TABLE]

We get thus

[TABLE]

According to Eqs. (11, 15, 12), we also obtain

[TABLE]

Taking the first derivative of $\phi(s,\lambda)$ from Eq. (125) with respect to $\lambda$ yields

[TABLE]

Since $\mu_{-}\mu_{+}=\sqrt{1-s}$, $\mu_{+}^{2}+\mu_{-}^{2}=2+s$, and $\mu_{+}+\mu_{-}=\sqrt{2+s+2\sqrt{1-s}}$, we get

[TABLE]

Expanding near $s\to 1$, one obtains

[TABLE]

from which

[TABLE]

in agreement with the general expansion (47).

Similarly, we compute $h_{2}(s)$ as the second derivative of $\phi(s,\lambda)$:

[TABLE]

As $s\to 1$, one gets

[TABLE]

from which

[TABLE]

in agreement with the general expansion (48).