MULTINOMIAL RANDOM COMBINATORIAL STRUCTURES AND

r-VERSIONS OF STIRLING, EULERIAN AND LAH NUMBERS

ALEXANDER IKSANOV, ZAKHAR KABLUCHKO, ALEXANDER MARYNYCH, VITALI WACHTEL

Abstract. We introduce multinomial and r-variants of several classic objects of combinatorial

probability, such as the random recursive and Hoppe trees, random set partitions and composi-
arXiv:2403.16448v1 [[Link]] 25 Mar 2024

tions, the Chinese restaurant process, Feller’s coupling, and some others. Just as various classic
combinatorial numbers – like Stirling, Eulerian and Lah numbers – emerge as essential ingredients
defining the distributions of the mentioned processes, the so-called r-versions of these numbers
appear in exact distributional formulas for the multinomial and r-counterparts. This approach
allows us to offer a concise probabilistic interpretation for various identities involving r-versions of
these combinatorial numbers, which were either unavailable or meaningful only for specific values
of the parameter r. We analyze the derived distributions for fixed-size structures and establish
distributional limit theorems as the size tends to infinity. Utilizing the aforementioned generalized
Stirling numbers of both kinds, we define and analyze (r, s)-Lah distributions, which have arisen in
the existing literature on combinatorial probability in various contexts.

Contents
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1. r-Stirling numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2. Three families of probability distributions related to the r-Stirling numbers . . . . . . . . 4
1.3. Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2. Multinomial structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1. Multinomial Stirling numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2. Multinomial Ewens permutation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3. Multinomial Chinese restaurant and the colored Feller coupling . . . . . . . . . . . . . . . . . . . . 9
2.4. Multinomial Hoppe forests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3. r-Ewens permutations, random r-recursive trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1. r-Ewens permutations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2. Random r-recursive trees and r-Hoppe trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4. Random r-partitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.1. Incomplete partitions derived from urn schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2. Random uniform r-partitions and random Gibbs r-partitions . . . . . . . . . . . . . . . . . . . . . . 22
5. Random r-compositions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.1. Couplings of random r-compositions and the Dirichlet multinomial distribution . . . . 26
5.2. Marginal distributions of r-compositions and their limits . . . . . . . . . . . . . . . . . . . . . . . . . . 29
6. The (r, s)-Lah distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
6.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2020 Mathematics Subject Classification. Primary: 60C05, 11B73; Secondary: 60E05, 60F05.
Key words and phrases. r-Stirling numbers, r-Lah numbers, central limit theorem, Chinese restaurant process,
Hoppe trees, random compositions, random permutations, random recursive trees, random set partitions.
1
 6.2. The r-Lah numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
6.3. Definition of the (r, s)-Lah distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6.4. Generating polynomial of the (r, s)-Lah distribution and a recursion for the
probability mass function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6.5. Stochastic representation in terms of random recursive trees . . . . . . . . . . . . . . . . . . . . . . . 34
6.6. Stochastic representation in terms of random compositions . . . . . . . . . . . . . . . . . . . . . . . . 36
6.7. Expectation and variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
7. Limit theorems for the (r, s)-Lah distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
Acknowledgement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

1. Introduction
1.1. r-Stirling numbers. The usual Stirling numbers of both kinds play a pervasive role in com-
binatorics and discrete probability, often appearing prominently in diverse counting problems as an
indispensable element. The list of quintessential examples begins with the enumeration of cycles in
permutations, records in random samples, and blocks in set partitions. It extends by enumerating
alignments and surjections, progressing to the counting of recursive trees with various profile re-
strictions. This list could have potentially been indefinitely supplemented. We refer to the excellent
treatises [10, 22, 46, 47, 50, 64] for comprehensive information on the classic Stirling numbers. A
thorough historic account and an extensive bibliography can be found in the books [46] and [47].
As detailed in the aforementioned references, there exist numerous generalizations of the Stirling
numbers. In the focus of the present paper are the so-called r-Stirling numbers which depend
on an additional parameter r ∈ R. These numbers have been intensively investigated in the
combinatorics literature. Given that the usual Stirling numbers play a central role in combinatorial
probability, it is natural to ask whether there exist natural r-deformations of the classic discrete
random structures (such as random permutations and partitions, or random recursive trees) whose
distributional properties are described by the r-Stirling numbers. One of the purposes of the present
paper is to construct these r-deformed structures and explore their properties.
Perhaps, the easiest way to define the r-Stirling numbers k r and nk r is by using the finite
n


differences calculus. Let C[r] be the vector space of polynomials in the variable r with coefficients
from C. Let D be the differentiation operator and ∆ the finite difference operator acting on C[r]
as follows:
Df (r) = f ′ (r), ∆f (r) = f (r + 1) − f (r).
Then, using the notation rn↑ = r(r + 1) · · · (r + n − 1) for the rising factorial,

Dk n↑ ∆k n
   
n n
:= r , := r , n ∈ N, k ∈ {0, 1, . . . , n}. (1.1)
k r k! k r k!

For example, n0 r = rn↑ and n0 r = rn . The usual Stirling numbers nk and nk are recovered by
     

evaluating the expressions at r = 0, that is,

     k       k 
n n D n↑ n n ∆ n
= = r , = = r .
k k 0 k! r=0 k k 0 k! r=0
2
In
n
general,
nthe parameter r may take arbitrary real (or complex) values. If we fix n and k, then
k r and k r are polynomials in r, as directly follows from (1.1). The above definitions can be
naturally extended to arbitrary n ∈ N0 := N ∪ {0} and k ∈ Z by putting
       
n n 0 0
:= := 0 for n ∈ N0 and k ∈
/ {0, . . . , n}, := := 1.
k r k r 0 r 0 r
The first systematic investigation of the r-Stirling numbers (also called the ‘weighted Stirling
numbers’ and ’non-central Stirling numbers’) is due to Carlitz [6, 7]. Some follow-up papers are [5,
12, 42, 58]. A clear exposition of properties of the r-Stirling numbers can be found in the paper
of Broder [5] (where the term ‘r-Stirling numbers’ appears for the first time) and in the books of
Mező [50, Chapter 8] and Charalambides [10, § 8.5], [11, Chapter 3]. The paper [54], whose notation
we follow, contains a nice overview of the identities satisfied by the r-Stirling numbers. We note  in
passing that our notation corresponds to Carlitz’s notation R(n, k, r) = nk r and R1 (n, k, r) = nk r


and what we denote by nk r and nk r corresponds to n+r

     n+r
k+r r and k+r r in [5, 50]. In particular,
the second definition in (1.1) can be read from Eq. (3.8) in [6].
While Carlitz’s approach is purely analytic and hinges on manipulations with generating func-
tions, Broder’s starting point is a combinatorial definition, see Eq. (1) and (2) in [5]. Adapting to
our notation and assuming that r is a non-negative integer, this is
Definition 1.1. For r ∈ N0 , n ∈ N and k ∈ {0, 1, . . . , n}, nk r is the number of permutations of
 

the set {1, 2, . . . , n + r} having k + r cycles such that the numbers 1, 2, . . . , r are in different cycles;
n

k r is the number of partitions of the set {1, 2, . . . , n + r} into k + r disjoint subsets such that the
numbers 1, 2, . . . , r are in distinct blocks.
While being combinatorially appealing Definition 1.1 has a drawback of being only valid for
integer r. In contrast, Carlitz’s method is free of this shortcoming. Unless stated otherwise, in all
subsequent considerations we stick to Carlitz’s viewpoint and do not assume that r is integer.
The properties of the r-Stirling numbers are parallel to those of the usual Stirling numbers. For
example, see Eq. (3.4) and (5.8) in [6] or Eq. (48) and (50) in [5],
n   n  
n↑
X n k n
X n
(x + r) = x , (x + r) = xk↓ , (1.2)
k r k r
k=0 k=0

where xk↓= x(x−1) · · · (x−k+1) denotes the falling factorial. The (vertical) exponential generating
functions of the r-Stirling numbers are given by
∞  
X n xn 1

k
= log (1 − x)−1 (1 − x)−r , (1.3)
k r n! k!
n=k
∞  
X n xn 1
= (ex − 1)k erx , (1.4)
k r n! k!
n=k
see Eq. (3.9) in [6] and Eq. (36) in [5]. It is known (see, e.g. p. 1661 in [54]) that
  n    n   
n X n j j−k X n j (n−j)↑
= r = r , (1.5)
k r j k j k
j=k j=k
  n    n   
n X n j n−j X n j (j−k)↓
= r = r . (1.6)
k r j k j k
j=k j=k
3
These formulas
n demonstrate
n once again that, for fixed n ∈ N0 , k ∈ {0, . . . , n}, the r-Stirling
numbers k r and k r are polynomials in r.

1.2. Three families of probability distributions related to the r-Stirling numbers. There
are various probability distributions related to usual Stirling numbers, see [11, 13, 53, 56, 57].
The most prominent example is the so-called Stirling distribution of the first kind (or a Stirling–
Karamata distribution) Stir1[n; θ] with parameters n ∈ N and θ > 0, which is a discrete probability
measure on {0, . . . , n} given by
  k
n θ
Stir1[n; θ]({k}) = , k ∈ {0, . . . , n}. (1.7)
k θn↑
If θ = 1, Stir1[n; 1] is the distribution of the number of cycles in a random permutation picked
uniformly at random from the symmetric group of n elements. More generally, Stir1[n; θ] is the
distribution of the number of cycles in a random Ewens permutation with parameter θ. This
distribution also pops up as the number of occupied tables in the Chinese restaurant process, see
Eq. (3.11) in [57], or the root degree in the random Hoppe tree [26]. For a review of the related
concept of the Ewens sampling formula and similar topics we refer to [15]. We will discuss these
connections later on in a larger generality.
Related to the Stirling numbers of the second kind is the Stirling–Sibuya distribution defined by
  k↓
n θ
StirSibuya[n; θ]({k}) = , k ∈ {0, . . . , n}, n ∈ N. (1.8)
k θn
Here, θ is either a positive integer or a real number larger than n − 1, so that the falling factorial is
nonnegative. If θ is a positive integer, StirSibuya[n; θ] is the distribution of the number of occupied
boxes after allocating n balls over θ urns equiprobably. The fact that both (1.7) and (1.8) do
define probability distributions follows from equalities (1.2) with r = 0. Yet another instance is the
Stirling distribution of the second kind Stir2[n; θ]. To define it, let Tn be the Touchard polynomial
given by
n  
X n k
Tn (z) = z , z ∈ C, n ∈ N.
k
k=0

Then, for a parameter θ > 0,

 
1 n k
Stir2[n; θ]({k}) = θ , k ∈ {0, . . . , n}, n ∈ N.
Tn (θ) k
This is the distribution of the number of blocks in a Gibbs (uniform, for θ = 1) random partition
of {1, . . . , n}. We refer to [11, 13, 37] and references therein for further properties and applications
of the aforementioned distributions.
From a purely analytic viewpoint one can introduce the r-versions of the above distributions by
simply using (1.2) and the r-Touchard polynomials.

Definition 1.2 ([53, 60, 61]). Let n ∈ N, θ > 0 and r ≥ 0. The r-Stirling distribution of the first
kind r-Stir1[n; θ] is a discrete probability measure on {0, . . . , n} given by

θk
 
n
r-Stir1[n; θ]({k}) = , k ∈ {0, . . . , n}. (1.9)
k r (θ + r)n↑
4
 The distribution defined by (1.9) appears in [53], see also [11, Example 2.1] and a survey [60].
Also, in [61] it pops up as the distribution of a random variable Xn+1 −1, where Xn+1 is the number
of records in a sample of size n + 1 of independent non-identically distributed random variables
(Nevzorov’s model). Our notation r-Stir1[n; θ] corresponds to Str1F(n, θ, r) in [53] and [61].
For r ≥ 0 and n ∈ N, let the r-Touchard polynomial, also called the r-Bell polynomial in [50,
§ 8.3] and [48], be defined by
n  
X n
Tn,r (z) := z k , z ∈ C.
k r
k=0

Definition 1.3. Let n ∈ N, θ > 0 and r ≥ 0. The r-Stirling distribution of the second kind
r-Stir2[n; θ] is a discrete probability measure on {0, . . . , n} given by
θk
 
n
r-Stir2[n; θ]({k}) = , k ∈ {0, . . . , n}. (1.10)
k r Tn,r (θ)
Definition 1.4 ([53, 60]). Assume that n ∈ N, r ≥ 0 and either θ ∈ N or θ > n − 1. The r-
Stirling–Sibuya distribution r-StirSibuya[n; θ] is a discrete probability measure on {0, . . . , n} given
by
θk↓
 
n
r-StirSibuya[n; θ]({k}) = , k ∈ {0, . . . , n}. (1.11)
k r (θ + r)n
The r-Stirling–Sibuya distribution r-StirSibuya[n; θ] appears in [53], see Section 5 therein, as
the distribution of a position at time n of a certain time-inhomogeneous random walk on N0 ; see
also [11, Example 2.6]. An urn model leading to the r-StirSibuya[n; θ] distribution is also discussed
in [53] and will be important in our considerations of random r-partitions in Section 4 below. Our
notation r-StirSibuya[n; θ] corresponds to Str2F(n, θ, r) in [53].
Distributions related to further generalizations of Stirling numbers are discussed in [11], [57,
Example 3.2.3] and [28].
In this work, we introduce a general concept of multinomial combinatorial structures which we
find interesting in their own. Multinomial structures are obtained by taking the classic random
structures (such as random permutations, partitions or compositions) and coloring their cycles or
blocks into d ∈ N colors. It will become clear that the r-distributions defined above (and some
other) occur in the special case d = 2.
The rest of the paper is organized as follows. In Section 2 we introduce and investigate the
multinomial versions of some classic combinatorial stochastic processes. In Sections 3 we show
how the r-distributions arise as the distributions of some functionals defined on the multinomial
structures of Section 2. Section 4 is devoted to the analysis of random incomplete partitions and
an accompanying notion of random r-partitions. A similar analysis of incomplete compositions and
random r-compositions is carried out in Section 5. In Section 6 we discuss yet another class of
combinatorial numbers, called r-Lah numbers, and introduce generalized Lah distributions. Spe-
cific occurrences of these distributions have already received some attention in the literature on
combinatorial probability across diverse contexts. Various limit theorems for these distributions
are discussed in the concluding Section 7.
1.3. Notation. Besides the already introduced notation, we also use the following abbreviations
and conventions. We put [n] := {1, 2, . . . , n}. For the asymptotic equivalence we use the symbol ≃.
Thus, an ≃ bn if, and only if, an /bn → 1 as n → ∞. We use the standard notation x+ = max(x, 0)
5
and x− = max(−x, 0) for x ∈ R. The cardinality of a finite set A is denoted by #A. For a formal
power series f (x) = ∞ n n
P
n=0 fn x , we stipulate that [x ]f (x) = fn . If X is a random element taking
values in some space S and µ is a probability distribution on S, we write X ∼ µ to denote that
X has distribution µ. If (µ[θ])θ∈P is a family of probability measures depending on a parameter θ
and M is a probability measure on P , we denote by µ[M] the mixture of (µ[θ])θ∈P . Thus,
Z
(µ[M])(·) := µ[θ](·)M(dθ). (1.12)
P

In addition to the already introduced probability measures related to the Stirling numbers, many
standard distributions show up throughout the text. We denote by Bin[n; p] the binomial distri-
bution with parameters n and p and by Bern[p] the Bernoulli distribution with parameter p, that
is, Bern[p]({1}) = p = 1 − Bern[p]({0}). Furthermore, Poi[λ] denotes the Poisson distribution
with parameter λ and Mult[n; (p1 , . . . , pd )] denotes the multinomial distribution with parameters
n and p1 , . . . , pd . Notation for other distributions will be introduced in the text just before the first
appearance. Keeping in mind Eq. (1.12) we ubiquitously use formulas like Bin[Poi[λ]; p] = Poi[λp]
involving the mixtures of distributions. With an abuse of notation, Bern[p] can denote either the
distribution or a random variable with this distribution (implicitly assumed to be defined on some
probability space). The same disclaimer also applies to the other distributions.

2. Multinomial structures
In this section we introduce r-analogues of some classic random structures. Examples include
random Ewens permutations, random recursive (and Hoppe) trees and the Chinese restaurant
process. The r-versions will be defined as special cases of more general, multinomial versions of
these structures.

 n 
2.1. Multinomial Stirling numbers. The multinomial Stirling cycle number k1 ,...,k d
is defined
as the number of permutations of [n] with k1 + · · · + kd cycles which are colored in d possible colors
1, . . . , d such that the number of cycles with color j is kj , for all j = 1, . . . , d. Clearly,
    
n n k1 + · · · + kd
= .
k1 , . . . , kd k1 + · · · + kd k1 , . . . , kd
 n
Similarly, the multinomial Stirling partition number k1 ,...,k d
is defined as the number of partitions
of [n] into k1 + · · · + kd nonempty blocks which are colored in d possible colors such that kj is the
number of blocks with color j, for all j = 1, . . . , d. Plainly,
    
n n k1 + · · · + kd
= .
k1 , . . . , k d k1 + · · · + kd k1 , . . . , kd

In both definitions,
 0the parameters n ∈ N0 and k1 , . . . , kd ∈ N0 satisfy k1 + · · · + kd ≤ n with the
convention that 0,...,0 = 1. Additionally, if n ∈ N0 and k1 , . . . , kd ∈ Z fail to satisfy the above
restrictions (meaning that one of the kj ’s is negative or their sum is larger than n), we define the
corresponding multinomial Stirling numbers to be 0. The properties of the multinomial Stirling
numbers are reviewed in [4, Section 4]. Here we only need the following result.
6
Proposition 2.1. With x1 , . . . , xd denoting variables,
 
X n
(x1 + · · · + xd ) = n↑
xk1 · · · xkdd , (2.1)
k1 , . . . , kd 1
k1 ,...,kd ∈N0
k1 +···+kd ≤n
 
n
xk11 ↓ · · · xkdd ↓ .
X
n
(x1 + · · · + xd ) = (2.2)
k1 , . . . , kd
k1 ,...,kd ∈N0
k1 +···+kd ≤n

Proof. To prove (2.1), we use the formula xn↑ = nk=0 nk xk with x = x1 + · · · + xd and then the
P  

multinomial theorem:
n   n    
X n X n X k
n↑
(x1 + · · · + xd ) = k
(x1 + · · · + xd ) = xk11 · · · xkdd
k k k1 , . . . , kd
k=0 k=0 k1 ,...,kd ∈N0
k1 +···+kd =k
    
X n k X n
= xk11 · · · xkdd = xk11 · · · xkdd .
k k1 , . . . , k d k1 , . . . , kd
k1 ,...,kd ∈N0 k1 ,...,kd ∈N0
k1 +···+kd ≤n k1 +···+kd ≤n
Pn n
Similarly, to prove (2.2), we use the formula xn = k=0 k xk↓ together with the analogue of the
multinomial theorem for falling/rising factorials:
n   n    
n n k
xk1 ↓ · · · xkdd ↓
X X X
n k↓
(x1 + · · · + xd ) = (x1 + · · · + xd ) =
k k k1 , . . . , kd 1
k=0 k=0 k1 ,...,kd ∈N0
k1 +···+kd =k
    
n k n
xk1 ↓ · · · xkdd ↓ = xk1 ↓ . . . xkdd ↓ .
X X
=
k k1 , . . . , k d 1 k1 , . . . , k d 1
k1 ,...,kd ∈N0 k1 ,...,kd ∈N0
k1 +···+kd ≤n k1 +···+kd ≤n

The proof is complete. □

2.2. Multinomial Ewens permutation. Let Sn be the set of all permutations of the set [n].
Let cycles(ρ) denote the set of cycles of a permutation ρ ∈ Sn and # cycles(ρ) denote the number
of cycles.
A random Ewens permutation with parameter θ > 0 is a random element ρ of Sn such that for
each deterministic permutation ρ ∈ Sn
θ# cycles(ρ)
P[ρ = ρ] = ,
θn↑
see Example 2.19 and Chapter 5 in [2], and also [57, Section 3.2]. We denote this distribution
by Ewens[n; θ], so that ρ ∼ Ewens[n; θ]. Note that the special case θ = 1 corresponds to a
random uniform permutation. It is important for what follows that the Ewens distribution can
also be defined for θ = 0. In this case it coincides with the uniform distribution on the subset of
Sn consisting of permutations with exactly one cycle. The total number of such permutations is
(n − 1)!.
The number of cycles of ρ follows the Stirling distribution of the first kind with parameter θ:
  k
n θ
P[# cycles(ρ) = k] = , k ∈ {0, . . . , n}
k θn↑
7
(in short, # cycles(ρ) ∼ Stir1[n; θ]). Indeed, nk is the number of permutations of [n] with k cycles,
 

and each such a permutation has probability θk /θn↑ under the Ewens distribution Ewens[n; θ].
Now we introduce the multinomial analogue of the Ewens distribution. A colored permutation
is a permutation of [n] in which each cycle has one of d possible colors. The number of colors
d ∈ N is fixed once and for all. Formally, a colored permutation is a pair (ρ, ϕ), where ρ ∈ Sn is a
permutation and ϕ : cycles(ρ) → {1, . . . , d} is a map assigning colors to cycles of ρ. A multinomial
random Ewens permutation is obtained by taking a Ewens permutation and coloring the cycles
independently of each other, with pj denoting the probability of assigning color j ∈ {1, . . . , d} to
a cycle. More precisely, we fix parameters θ > 0 and p1 , . . . , pd ≥ 0 such that p1 + · · · + pd = 1,
define θj := θpj and adopt the following

Definition 2.2. A random colored permutation (ρ, ϕ) has a multinomial Ewens distribution
MultiEwens[n; θ; (p1 , . . . , pd )] if
# cycles1 (ρ,ϕ) # cyclesd (ρ,ϕ)
θ# cycles(ρ) # cycles1 (ρ,ϕ) # cyclesd (ρ,ϕ) θ . . . θd
P[(ρ, ϕ) = (ρ, ϕ)] = n↑
p1 . . . pd = 1 (2.3)
θ θn↑
for each deterministic colored permutation (ρ, ϕ). Here, cyclesj (ρ, ϕ) denotes the set of cycles of
color j in the colored permutation (ρ, ϕ).

From the definition of multinomial Stirling numbers it follows that the cycle counts of such (ρ, ϕ)
are distributed as follows:
  k1 +···+kd
n θ
P[# cycles1 (ρ, ϕ) = k1 , . . . , # cyclesd (ρ, ϕ) = kd ] = pk11 · · · pkdd
k1 , . . . , kd θn↑
for all k1 , . . . , kd ∈ N0 subject to k1 + · · · + kd ≤ n. This motivates the following

Definition 2.3. A random vector (K1 , . . . , Kd ) with values in Nd0 has a multinomial Stirling dis-
tribution of the first kind MultiStir1[n; θ; (p1 , . . . , pd )] if

θ1 · · · θdkd
  k1 +···+kd   k1
n θ k1 kd n
P[K1 = k1 , . . . , Kd = kd ] = p1 · · · pd =
k1 , . . . , k d θn↑ k1 , . . . , kd θn↑
for all k1 , . . . , kd ∈ N0 satisfying k1 + · · · + kd ≤ n, where θj = θpj for j ∈ {1, . . . , d}.

Thus, the cycle counts of (ρ, ϕ) ∼ MultiEwens[n; θ; (p1 , . . . , pd )] satisfy

(# cycles1 (ρ, ϕ), . . . , # cyclesd (ρ, ϕ)) ∼ MultiStir1[n; θ; (p1 , . . . , pd )].

Proposition 2.4. The multivariate generating function of the random vector (K1 , . . . , Kd ) that
follows MultiStir1[n; θ; (p1 , . . . , pd )] distribution is given by

Kd (θ1 t1 + · · · + θd td )n↑
E[tK
1 · · · td ] =
1
, (2.4)
(θ1 + · · · + θd )n↑
where θj = θpj for i ∈ {1, . . . , d}. Moreover, if I1 , . . . , In are independent random vectors in Nd0
with
θj ℓ−1
P[Iℓ = ej ] = , j ∈ {1, . . . , d}, P[Iℓ = 0] = , (2.5)
θ+ℓ−1 θ+ℓ−1
where e1 , . . . , ed is the standard basis of Rd and 0 is the zero vector in Rd , then (K1 , . . . , Kd ) has
the same distribution as I1 + · · · + In . In particular, Kj ∼ (θ − θj )-Stir1[n; θj ].
8
Proof. By definition of the multinomial Stirling distribution,
(θ1 t1 + · · · + θd td )n↑
 
K1 Kd 1 X n
E[t1 · · · td ] = n↑ (θ1 t1 )k1 . . . (θd td )kd = ,
θ k1 , . . . , k d θn↑
k1 ,...,kd ∈N0
k1 +···+kd ≤n

where the second equality follows from (2.1). It remains to note that the generating function of
the random vector I1 + · · · + In is the same. To identify the distribution of Kj , observe that (2.4)
implies
n  
Kj (θj tj + θ − θj )n↑ (1.2) X n θjk tkj
E[tj ] = = ,
θn↑ k θ−θj (θj + θ − θj )n↑
k=0
which proves that Kj ∼ (θ − θj )-Stir1[n; θj ], see Definition 1.2. □
The following aggregation property follows directly from the definition:

Proposition 2.5. If (K1 , . . . , Kd ) ∼ MultiStir1[n; θ; (p1 , p2 , . . . , pd )], then

(K1 + K2 , K3 , . . . , Kd ) ∼ MultiStir1[n; θ; (p1 + p2 , p3 , . . . , pd )].

Recalling our convention (1.12) concerning the mixtures of distributions we write

MultiStir1[n; θ; (p1 , . . . , pd )] = Mult[Stir1[n; θ]; (p1 , . . . , pd )]. (2.6)
This provides another way to prove Proposition 2.4.

2.3. Multinomial Chinese restaurant and the colored Feller coupling. Fix parameters
θ1 , . . . , θd ≥ 0 and suppose that θ := θ1 + · · · + θd > 0. Customers labeled by 1, 2, . . . , n arrive in
discrete time and take seats at round tables. Each table is colored in one of d colors. Each customer
ℓ, after arrival, either creates a new table of color j ∈ {1, . . . , d} with probability θj /(θ + ℓ − 1) or
takes a seat at some already existing table next to customer i (in the counterclockwise direction)
with probability 1/(θ + ℓ − 1), for each customer i = 1, . . . , ℓ − 1. After arrival of n customers,
this process generates a random permutation of [n], if we regard tables (read counterclockwise) as
cycles.

Proposition 2.6. The colored permutation generated by the multinomial Chinese restaurant process
has the distribution
MultiEwens[n; θ; (θ1 /θ, . . . , θd /θ)].
In particular, if Cn;j denotes the number of tables (cycles) of color j ∈ {1, . . . , d}, then
(Cn;1 , . . . , Cn;d ) ∼ MultiStir1[n; θ; (θ1 /θ, . . . , θd /θ)].

First proof. Consider the multinomial Chinese restaurant process and remove the colors of the
tables. This is the usual Chinese restaurant process with parameter θ. It is known, see Example
2.19 in [2], that the process generates a random permutation which is distributed according to
Ewens[n; θ]. Now, recall that in the colored version of the process each new table receives color j
with probability θj /θ. By definition, this gives the multinomial Ewens distribution. □
Second proof. Another approach hinges on a colored version of the classic Feller coupling for Ewens’
permutations, see pp. 95-96 in [2]. Let (Di )i∈N , be independent random variables such that Di takes
values in the set
{close1 , . . . , closed } ∪ {2, 3, . . . , i}
9
(the meaning of ‘closej ’ will be explained a bit later) and has the distribution
θj 1
P[Di = closej ] = , j = 1, . . . , d, P[Di = k] = , k = 2, . . . , i, i ∈ N.
θ+i−1 θ+i−1
These variables are then used to generate a colored random permutation in the ordered cycle
notation as follows. Start the first cycle with ‘(1’ and make a (n + d − 1)-choice using the variable
Dn . If Dn = k, then k is added to the first cycle transforming ‘(1’ to ‘(1 k’. If Dn = closej , then
one closes the first cycle, paints it using the color j and opens the next cycle. Continuing in this
way and using Dn−1 , . . . , D1 produces a colored random permutation, say (ρ, ϕ), in the ordered
cycle notation. Thus, the value Di = closej corresponds to closure of the current cycle, painting it
using the color j and starting a new cycle with the smallest unused integer. The total number of
cycles in (ρ, ϕ) is equal to
Xn
Cn = 1[Di ∈ {close1 , . . . , closed }]
i=1
and, since P[Di ∈ {close1 , . . . , closed }] = θ/(θ + i − 1), follows the Stirling distribution of the first
kind Stir1[n; θ]. Furthermore, the probability of seeing any fixed (uncolored) permutation ρ ∈ Sn ,
is equal to
θ#cycles(ρ)
P[ρ = ρ] = .
θn↑
Thus, erasing the colors in (ρ, ϕ) produces Ewens’ permutation. The number of cycles of color j is
n
X
Cn;j = 1[Di = closej ], j = 1, . . . , d.
i=1
Using that, independently of i = 1, . . . , n,
θj
P[Di = closej |Di ∈ {close1 , . . . , closed }] = , j = 1, . . . , d,
θ
we conclude that (ρ, ϕ) ∼ MultiEwens[n; θ; (θ1 /θ, . . . , θd /θ)] and (Cn;1 , . . . , Cn;d ) has the multino-
mial Stirling distribution MultiStir1[n; θ; (θ1 /θ, . . . , θd /θ)]. □
Since the cycles are painted independently, the following is a simple consequence of the marking
theorem for Poisson processes in conjunction with Theorem 5.1 in [2].

Proposition 2.7. Let (ρ, ϕ) be the multinomial Ewens permutation and let Cn;j;k be the number
of cycles of color j ∈ {1, . . . , d} and length k ∈ N. Then,
d
(Cn;j;k : j ∈ {1, . . . , d}, k ∈ N) −→ (Poi[θj /k] : j ∈ {1, . . . , d}, k ∈ N)
n→∞

in the product topology on Nd×N

0 , where the limiting Poisson variables are mutually independent.

2.4. Multinomial Hoppe forests. The classic sequential construction of a random recursive
tree RRT[n] with n nodes 1, . . . , n and the root labeled by 0 proceeds by recursively attaching
node ℓ ∈ {1, . . . , n} to one of the existing nodes 0, 1, . . . , ℓ − 1 picked uniformly at random and
independently of the previous choices; see, e.g., [62]. The Hoppe tree Hoppe[n; θ] is defined similarly,
see [43], but now the probability of attaching the node ℓ to the root is θ/(θ + ℓ − 1), while the
probability of attaching ℓ to any other node 1, . . . , ℓ − 1 is 1/(θ + ℓ − 1). Here, θ > 0 is a parameter
representing the ‘weight’ of the root, while the ‘weight’ of any other node is 1. In the special case
θ = 1 one recovers the random recursive tree RRT[n]. The name ‘Hoppe tree’ originates from the
10
corresponding Hoppe’s urn model introduced in [25, 26, 27] in the context of the Ewens sampling
formula.
We now define a multinomial generalization of Hoppe trees, called multinomial Hoppe forest and
denoted by MultiHoppe[n; θ; (θ1 /θ, . . . , θd /θ)]. Fix parameters θ1 , . . . , θd ≥ 0 and suppose that
θ := θ1 + · · · + θd > 0. Start at time 0 with d roots, labeled root1 , . . . , rootd , which have weights
θ1 , . . . , θd . Nodes, labeled by 1, 2, . . . , n, arrive in discrete time and are added as children to the
roots (with probabilities proportional to the weights θ1 , . . . , θd ) or to the nodes already present
in the forest (with probability proportional to 1 for each node). More precisely, suppose that the
nodes 1, . . . , ℓ − 1 have already arrived, where ℓ ∈ {1, . . . , n}. Then, the next node, with label ℓ,
arrives at time ℓ and is attached to rootj with probability θj /(θ + ℓ − 1), for each j ∈ {1, . . . , d}, or
to any of the already existing nodes 1, . . . , ℓ − 1, with probability 1/(θ + ℓ − 1) for each node. The
usual Hoppe trees are recovered by taking d = 1.
Recall that a random vector (X1 , . . . , Xd ) has the Dirichlet multinomial distribution, denoted
hereafter by DirMult[n; θ1 , . . . , θd ], if

P[(X1 , . . . , Xd ) = (k1 , . . . , kd )]
d d
Γ(n + 1)Γ(θ1 + · · · + θd ) Y Γ(ki + θi ) X
= , ki ≥ 0, ki = n, (2.7)
Γ(n + θ1 + · · · + θd ) Γ(θi )Γ(ki + 1)
i=1 i=1

see Eq. (35.151) in [35] with the convention that P[Xi = 0] = 1 if θi = 0, for some i ∈ {1, . . . , d}.
Let Beta[α; β] denote the beta distribution with parameters α, β > 0 and having the density
x 7→ xα−1 (1 − x)β−1 / B(α, β) for x ∈ (0, 1). If α = 0 (respectively β = 0), we stipulate that
Beta[α; β] is the degenerate distribution at 0 (respectively, 1). The following result is standard,
see [34, Chapter 4.5.1].

Proposition 2.8. The joint distribution of the sizes of the connected components of root1 , . . . , rootd
in MultiHoppe[n; θ; (θ1 /θ, . . . , θd /θ)] (not counting the roots themselves) is DirMult[n; θ1 , . . . , θd ].
In particular, the size of the connected component of rootj (not counting rootj itself ) has the so-
called beta-binomial distribution Bin[n; Beta[θj ; θ − θj ]], for all j = 1, . . . , d.

Proposition 2.9. For each j ∈ {1, . . . , d}, let Dn;j be the degree of rootj in the random forest
MultiHoppe[n; θ; (θ1 /θ, . . . , θd /θ)]. Then,
(Dn;1 , . . . , Dn;d ) ∼ MultiStir1[n; θ; (θ1 /θ, . . . , θd /θ)].

Proof. For ℓ ∈ {1, . . . , n} let Iℓ be a random vector in Nd0 with distribution (2.5) and assume that
I1 , . . . , In are mutually independent. We interpret the event {Iℓ = ej } as {the node ℓ was attached
to rootj } and the event {Iℓ = 0} as {ℓ was attached to one of the nodes 1, . . . , ℓ − 1}. Then,
(Dn;1 , . . . , Dn;d ) has the same distribution as I1 + · · · + In . It remains to apply Proposition 2.4. □

2.4.1. Expected profiles of multinomial Hoppe forests.

  The expected number of nodes in an RRT[n]
1 n+1
having distance k to the root is known to be n! k+1 , for all k ∈ {1, . . . , n}; see [18, Lemma 6.16]
or [62, Theorem 1] (actually both references attribute this formula to [17]). The next theorem gen-
eralizes this result by identifying the expected profile of the multinomial Hoppe trees. Interestingly,
even for d = 1 and arbitrary θ > 0, the next formula contains the θ-Stirling numbers.

Theorem 2.10. Consider MultiHoppe[n; θ; (θ1 /θ, . . . , θd /θ)] and let Ln;j (k) be the number of
nodes in the connected component of the root j ∈ {1, . . . , d} which have distance k to the root
11
j. Then,  
θj n
E[Ln;j (k)] = n↑ , k ∈ {1, . . . , n}, j ∈ {1, . . . , d}. (2.8)
θ k θ
Proof. Denote the left-hand side of (2.8) by aj (n, k) := E[Ln;j (k)]. For n ∈ N0 , let Fn be the
σ-algebra generated by the multinomial Hoppe forest, which starts with d roots, after adding n
new nodes. Then aj (0, k) = 0 for k ∈ N and, for k ≥ 2,
 
Ln−1;j (k − 1)
aj (n, k) = E[Ln;j (k)] = E[E[Ln;j (k)]|Fn−1 ] = E Ln−1;j (k) +
n+θ−1
aj (n − 1, k − 1)
= aj (n − 1, k) + , n ∈ N, j = 1, . . . , d, (2.9)
n+θ−1
where 1/(n + θ − 1) is the probability that the node n joins a particular node at the level k − 1.
For k = 1, the recursion is slightly different accounting to a special weight of the root j:
θj
aj (n, 1) = aj (n − 1, 1) + , n ∈ N, j = 1, . . . , d, (2.10)
n+θ−1
where θj /(n + θ − 1) is the probability of joining the root j. Introduce the generating functions
∞
X
ϕn;j (z) := θj + aj (n, k)z k , n ∈ N0 , j = 1, . . . , d.
k=1

Then (2.9) and (2.10) imply

 
z
ϕ0;j (z) = θj , ϕn;j (z) = 1+ ϕn−1;j (z), n ∈ N.
n+θ−1
Thus,
n−1 n  
Y (θ + z)n↑

z θj X n k
ϕn;j (z) = θj 1+ = θj = z , n ∈ N0 ,
k+θ θn↑ θn↑ k θ
k=0 k=0
and (2.8) follows. □

Remark 2.11. Observe that the expected number of nodes in the j-th tree (not counting the root)
is
n n  
X θj X n (1.2) θj
E[Ln;j (k)] = n↑ = n↑ ((1 + θ)n↑ − θn↑ ) = nθj /θ,
θ k θ θ
k=1 k=1
n
as it should be. We have used here that 0 θ = θn↑ .

2.4.2. Leaves in multinomial Hoppe forests and generalized Eulerian numbers. Let Tn be the number
of leaves in a random recursive tree RRT[n] with one root and n other nodes. According to a result
of Najock and Heyde [52] (alternatively, see Eq. (6.15) in [18]):
 
1 n
P[Tn = k] = , k = 1, . . . , n, n ∈ N, (2.11)
n! k − 1
where the second factor is the standard Eulerian number, see [22, Chapter 6.2] and [55]. For our
purposes the following recursive definition of the Eulerian number serves best:
     
n n−1 n−1
= (n − k) + (k + 1) , 0 ≤ k < n, n ∈ N (2.12)
k k−1 k
12
with the boundary conditions n0 = 1, n ∈ N0 , nn = 0, n ∈ N. In the following we generalize (2.11)
first to Hoppe trees and then to Hoppe forests. Put A(r, s) := r+s+1
s and note that (2.12) implies
A(r, s) = (r + 1)A(r, s − 1) + (s + 1)A(r − 1, s), r, s ∈ N0 , r + s > 0,
with A(0, 0) = 1 and A(r, s) = 0 if r = −1 or s = −1. This formula suggests the following
generalization of Eulerian numbers, proposed by Carlitz and Scoville [8]. Let α, β be arbitrary real
parameters. Put
A(r, s|α, β) = (r + β)A(r, s − 1|α, β) + (s + α)A(r − 1, s|α, β), r, s ∈ N0 , r + s > 0, (2.13)
and A(0, 0|α, β) = 1 and A(r, s|α, β) = 0 if r = −1 or s = −1. Note that A(r, s) = A(r, s|1, 1).
According to Eq. (1.8) in [8],
∞
X ur v s
A(r, s|α, β) = (1 + uF (u, v))α (1 + vF (u, v))β ,
(r + s)!
r,s=0

where
eu − ev
F (u, v) = .
uev − veu
It turns out that the generalized Eulerian numbers defined by (2.13) are an indispensable ingredient
of the version of (2.11) for Hoppe trees.

Theorem 2.12 (Leaves in Hoppe trees). Let Tn be the number of leaves in a random tree Hoppe[n; θ]
with θ > 0. Then,
A(n − k, k|0, θ)
P[Tn = k] = , k = 1, . . . , n, n ∈ N. (2.14)
θn↑
Proof. Denote by pn (k) := P[Tn = k] the left-hand side of (2.14) and let Fn denote the σ-algebra
generated by the Hoppe tree upon adding n nodes. Note that the root cannot be a leaf if n ≥ 1
and it is also convenient not to regard it as a leaf if n = 0, thus we put p0 (0) = 1. Then,

k n−k+θ
pn (k) = E[P[Tn = k|Fn−1 ]] = pn−1 (k)+ pn−1 (k−1), k = 0, . . . , n, n ∈ N,
n+θ−1 n+θ−1
where we stipulate that pn (k) = 0, for k > n, and pn (−1) = 0. Put Qn;ℓ := θ(n+ℓ)↑ pn+ℓ (ℓ) and note
that the above recursion transforms into
Qn;ℓ = ℓQn−1;ℓ + (n + θ)Qn;ℓ−1 , n ∈ N0 , ℓ ∈ N0 , n+ℓ>0
with the initial values
Q0,0 = 1, Q−1;ℓ = Qn;−1 = 0, n, ℓ ∈ N0 .
Thus, Qn;ℓ = A(n, ℓ|0, θ) and (2.14) follows. □

Example 2.13. Note that the RRT result (2.11) is a particular case when θ = 1, since A(n −
n r+s
k, k|0, 1) = k−1 which, in turn, is a consequence of A(r, s|0, 1) = s−1 .

The generalized Eulerian distribution with probability mass function k 7→ A(k, n − k|α, β)/(α +
β)n↑ for k ∈ {0, . . . , n} is investigated in [9, 30, 31] in the context of Morisita’s work on mutual
repulsive behavior of ant lions [51]. Among other results, Charalambides [9] computes the factorial
moments of this distribution and proves a central limit theorem.
13
Remark 2.14. The r-analogues of Eulerian numbers appearing in [50, § 8.5–8.7] and [49] and denoted
there by nk r satisfy the same recurrence relation, namely, nk r = (n − k + r) n−1 n−1
k−1 r + (k + 1) k r
(see relation (8.30) on p. 218 in [50]) as the numbers A(n − k − 1, k + 1|0, r + 1) (see (2.13)).
However, the boundary conditions imposed on these arrays are different; see [50, p. 216]. Related
generalizations of Eulerian numbers can be found in [32] and [45].

Theorem 2.15 (Leaves in Hoppe forests). Let Tn;j be the number of leaves in the connected compo-
nent of the root j ∈ {1, . . . , d} in the multinomial Hoppe forest MultiHoppe[n; θ; (θ1 /θ, . . . , θd /θ)]
(if the root j has no descendants, we stipulate that Tn,j = 0). Then,
n  
1 X n
P[Tn;j = k] = n↑ (θ − θj )(n−m)↑ A(m − k, k|0, θj ).
θ m
m=k

Proof. The number of nodes in the connected component of rootj (not counting the root itself) is
Sn;j ∼ Bin[n; Beta[θj ; θ − θj ]], see Proposition 2.8. Conditionally on Sn;j = m, the component
of rootj is distributed as Hoppe[m; θj ]. By the total probability formula and Theorem 2.12, the
number of leaves in the connected component of rootj satisfies
n
X A(m − k, k|0, θj )
P[Tn;j = k] = P[Sn;j = m] ·
m=k θjm↑
n  
X n Γ(m + θj )Γ(n − m + θ − θj )Γ(θ) Γ(θj )
= · A(m − k, k|0, θj )
m Γ(n + θ)Γ(θ − θj )Γ(θj ) Γ(m + θj )
m=k
n  
1 X n
= n↑ (θ − θj )(n−m)↑ A(m − k, k|0, θj ),
θ m
m=k

which proves the claim. □

Corollary 2.16 (Subtrees of the ℓ-th node in a Hoppe tree). Consider a Hoppe[n; θ]-tree and let
(ℓ)
ℓ ∈ {1, . . . , n} be a node. If Tn denotes the number of leaves in the subtree rooted at the node ℓ,
then
n−ℓ    
1 X n−ℓ m
P[Tn(ℓ) = k] = (ℓ + θ − 1)(n−ℓ−m)↑
k = 1, . . . , n − ℓ
(θ + ℓ)(n−ℓ)↑ m=k m k−1

(ℓ)
and P[Tn = 0] = (ℓ − 1 + θ)/(n − 1 + θ), which is the probability that the subtree only consists of ℓ.

Proof. Indeed, after the node ℓ has arrived, we can consider it as root1 with weight θ1 = 1 and
combine the nodes 1, . . . , ℓ − 1 and the old root to form root2 with weight θ2 = ℓ − 1 + θ. The
remaining nodes ℓ + 1, . . . , n behave as in the corresponding multinomial Hoppe tree with d = 2.
m
Applying Theorem 2.15 and recalling that A(m − k, k|0, 1) = k−1 gives the claim. □

In the case of RRT (that is, for θ = 1) our formula is equivalent to the formula of Mahmoud and
Smythe [44, p. 410]. Note, however, that they stipulate a tree consisting of a single root to have
one leaf, whereas we declare it to have no leaves. Summarizing, their formula gives the distribution
(ℓ)
of max(Tn−1 , 1).
14
 3. r-Ewens permutations, random r-recursive trees
Now we are going to define the r-analogues of the random structures introduced above. The
general procedure is as follows. Fix parameters τ = θ1 ≥ 0 and r = θ2 ≥ 0 such that θ = τ + r > 0.
Consider some multinomial random structure with d = 2 colors having probabilities p1 := τ /(τ + r)
and p2 := r/(τ + r) and discard all elements having color 2. The remaining elements of color 1 form
a random structure on some subset B ⊆ [n] (which is also random and, moreover, can be empty
with positive probability). This random structure, defined on a random subset, is the r-analogue
we are searching for.

3.1. r-Ewens permutations. The r-Ewens distribution r-Ewens[n; τ ] is informally defined as

follows. Consider a random permutation following the Ewens[n; τ + r] distribution. Independently
remove each cycle of this permutation with probability r/(τ + r). The remaining cycles define a
random permutation of some subset B of [n], which we call the r-Ewens permutation and whose
distribution is r-Ewens[n; τ ]. We now give a precise definition.

Definition 3.1. An incomplete permutation on [n] is a pair (ρ, B), where B ⊆ [n] is a subset of
[n], which might be empty, and ρ : B → B is a bijection. The set B is called the domain of ρ and
its elements is colored white. The complement of B is denoted by B0 := [n]\B and its elements
are colored red.

The set of incomplete permutations of [n] is denoted by Sinc

n = {(ρ, B) : B ⊆ [n], ρ ∈ SB }, where
SB is the symmetric group of B.

Definition 3.2. An r-Ewens permutation with parameters τ ≥ 0 and r ≥ 0 such that τ + r > 0 is
a random variable (ρ, B) taking values in Sinc
n with the distribution

τ # cycles(ρ)
P[(ρ, B) = (ρ, B)] = r(n−#B)↑ · , (ρ, B) ∈ Sinc
n , (3.1)
(τ + r)n↑
which is called the r-Ewens[n; τ ] distribution. If τ = 1, then we say that (ρ, B) is an r-uniform
random permutation.

Proposition 3.3. Consider a multinomial Chinese restaurant process with d = 2 colors having
weights θ1 = τ and θ2 = r. Let (ρ, B) be an incomplete permutation whose cycles are the tables of
color 1 and whose domain B is the set of all customers sitting at tables of color 1. Let B0 := [n]\B
be the set of all customers sitting at tables of color 2. Then, (ρ, B) ∼ r-Ewens[n; τ ].

Proof. Fix (ρ, B) ∈ Sincn and let Ωρ,B be the set of colored permutations of [n] (with two colors)
such that the elements belonging to cycles with color 2 are precisely the elements of the set B0
and cycles of color 1 form permutation ρ of B. According to Proposition 2.6 we need to check that
the sum of the right-hand sides of (2.3) over the set Ωρ,B is equal to the right-hand side of (3.1).
Indeed,
′ ′ ′
X τ # cycles1 (ρ ,ϕ) r# cycles2 (ρ ,ϕ) τ # cycles(ρ) (n−#B)↑ X r# cycles2 (ρ ,ϕ)
= r .
(t + r)n↑ (t + r)n↑ r(n−#B)↑
(ρ′ ,ϕ)∈Ωρ,B (ρ′ ,ϕ)∈Ωρ,B

The sum is equal to 1, since the summands form a probability distribution Ewens[n − #B; r] and
the summation is taken over all possible cyclic structures on [n] \ B. □
15
Proposition 3.4. Let An,r (b0 , b1 , . . . , bk ) be the event that (ρ, B) ∼ r-Ewens[n; τ ] has a red set of
size b0 ∈ N0 and k white cycles of sizes b1 , . . . , bk ∈ N, where b0 + b1 + · · · + bk = n. Then,
n! τk rb0 ↑ 1
P[An,r (b0 , b1 , . . . , bk )] = n↑
· · . (3.2)
k! (τ + r) b0 ! b1 · · · bk
n

1
Proof. Equality (3.2) follows immediately from (3.1) upon noticing that b0 ,...,b k k!
counts the num-
ber of partitions of [n] into blocks B0 , B1 , . . . , Bk of sizes b0 ∈ N0 , b1 , . . . , bk ∈ N, respectively,
whereas (bj − 1)! is the number of ways to organize a cycle on bj elements, j = 1, . . . , k. □
Proposition 3.5. If (ρ, B) is an r-Ewens[n; τ ]-distributed incomplete permutation, then the num-
ber of cycles in ρ has the r-Stir1[n; τ ] distribution, that is,
τk
 
n
P[# cycles(ρ) = k] = , k ∈ {0, . . . , n}. (3.3)
k r (τ + r)n↑
Before proceeding with a proof of Proposition 3.5 we note that upon multiplying Taylor expan-
sions (1.3) and (1.4), one obtains
rb0 ↑ r b0
   
n n! X 1 n n! X 1
= · , = · , (3.4)
k r k! b1 · · · bk b0 ! k r k! b1 ! · · · bk ! b0 !
(b0 ,b1 ,...,bk ) (b0 ,b1 ,...,bk )

where the sums are taken over all (k + 1)-tuples (b0 , b1 , . . . , bk ) such that b0 ∈ N0 , b1 , . . . , bk ∈ N
and b0 + b1 + · · · + bk = n.
We give three elementary proofs of Proposition 3.5 with the aim of clarifying the probabilistic
meaning of the first equality in (3.4) (the second will be explained later on, see Proposition 4.5)
and also of the formulas
  n    n   
n X n j j−k X n j (n−j)↑
= r = r , (3.5)
k r j k j k
j=k j=k

see Eq. (1.5) in the introduction.

First proof. Formula (3.3) for the distribution of the number of cycles follows from the first equality
in (3.4) upon summation of the right-hand sides of (3.2) over all b0 ∈ N0 and b1 , . . . , bk ∈ N satisfying
b0 + · · · + bk = n. □
Second proof. We use the first equality in (3.5). Since the total number of cycles of both colors
has the Stir1[n; τ + r] distribution and the probability for each particular cycle to be of color 1 is
τ /(τ + r) independently of the other cycles,
P[# cycles(ρ) = k] = Bin[Stir1[n; τ + r]; τ /(τ + r)]({k}), k ∈ {0, . . . , n}.
Thus, we need to check that
r-Stir1[n; τ ] = Bin[Stir1[n; τ + r]; τ /(τ + r)]. (3.6)
This equality follows from
n
X
Bin[Stir1[n; τ + r]; τ /(τ + r)]({k}) = Stir1[n; τ + r]({j})Bin[j; τ /(τ + r)]({k})
j=k
n  
n (τ + r)j j τ k rj−k (3.5) n τk
X    
= = = r-Stir1[n; τ ]({k})
j (τ + r)n↑ k (τ + r)j k r (τ + r)n↑
j=k
16
for each k ∈ {0, 1, . . . , n}. □
Third proof. In this proof we exploit the second equality in (3.5). Using definition (3.1) of the
r-Ewens distribution we obtain
n n (n−j)↑ τ k n
  
X (3.1) X r j
P[# cycles(ρ) = k] = P[# cycles(ρ) = k, #B = j] = n↑
, k ∈ {0, . . . , n},
(τ + r) j k
j=k j=k

where nj counts the number of ways to choose j elements of [n] and kj counts the number of ways

 

to organize k cycles on the chosen elements. By the second equality in (3.5) the right-hand side of
the last centered formula is equal to r-Stir1[n; τ ]({k}). □
Similarly to the classic scenario, the number of cycles in the r-Ewens permutation (ρ, B) can be
represented as the sum of independent indicators. This follows from the factorization
n   n 
(r + τ s)n↑

# cycles(ρ) 1 X n k
Y r+j−1 τs
E[s ]= (τ s) = = + .
(r + τ )n↑ k r (r + τ )n↑ r+τ +j−1 r+τ +j−1
k=0 j=1

This representation immediately yields a central limit theorem for # cycles(ρ) and, therefore, for
the r-Stirling distribution of the first kind. Interestingly, the limit theorem does not depend on r
at all!
Corollary 3.6. Assume that τ, r ≥ 0, τ +r > 0 are fixed and let (ρ, B) be the r-Ewens permutation.
Then,
# cycles(ρ) − τ log n d
√ −→ N [0; 1] .
τ log n n→∞

For the ease of reference we collect various representations for the r-Stirling distribution of the
first kind derived above in a single
Proposition 3.7. Fix r ≥ 0, τ > 0. The r-Stirling distribution of the first kind r-Stir1[n; τ ] admits
the following representations:
Pn
(i) j=1 Bern[τ /(τ + r + j − 1)] ∼ r-Stir1[n; τ ], where the random variables on the left-hand
side are independent;
(ii) r-Stir1[n; τ ] = Bin[Stir1[n; τ + r]; τ /(τ + r)];
(iii) r-Stir1[n; τ ] = Stir1[Bin[n; Beta[τ ; r]]; τ ], where Beta[τ ; r] is the beta-distribution.
Proof. Part (i) can be verified by comparing generating functions. Alternatively, it is just the last
claim of Proposition 2.4 with d = 2, j = 1, θ1 = τ , θ2 = r. Part (ii) is formula (3.6). For the proof
of part (iii) note that, for k ∈ {0, 1, . . . , n},
n
X
Stir1[Bin[n; Beta[τ ; r]]; τ ]({k}) = Bin[n; Beta[τ ; r]]({m})Stir1[m; τ ]({k})
m=0
n 
n Γ(m + τ )Γ(n − m + r) Γ(τ + r) m τ k
X   
=
m Γ(n + τ + r) Γ(τ )Γ(r) k τ m↑
m=0
n   
Γ(τ + r) k
X n m (n−m)↑ Γ(m + τ )
= τ r
Γ(n + τ + r) m k Γ(τ )τ m↑
m=0
τk
 
(1.5) n
= n↑
= r-Stir1[n; τ ]({k}).
(τ + r) k r
17
The proof is complete. □
Remark 3.8. The recent article [29] investigates properties of the Bernoulli trials with unequal
harmonic-like success probabilities as in part (i) of Proposition 3.7.
We close the discussion of r-Ewens permutations by noting a combinatorial construction of the
r-uniform permutations which works for integer r ∈ N0 and establishes connections to Broder’s
Definition 1.1. Call the elements 1, 2, . . . , n of the set [n + r] white and elements n + 1, . . . , n + r
red. Let Sn,r be the set of all permutations of [n + r] such that red elements are in different cycles
and let dperm : Sn,r 7→ Sinc
n be the mapping defined as follows. For σ ∈ Sn,r , the image dperm (σ) is
an incomplete permutation (ρ, B) ∈ Sinc n such that the red set B0 = [n] \ B is obtained by merging
r cycles of σ containing red elements, removing these red elements from the union, and keeping
other cycles unchanged.
Proposition 3.9. Fix n ∈ N, r ∈ N0 . Let σ bn,r be a permutation of [n + r] picked uniformly at
random from Sn,r . Then dperm (b
σn,r ) is the r-uniform random permutation (r-Ewens permutation
with τ = 1).
Proof. Each σ ∈ Sn,r can be constructed from an incomplete permutation (σ1 , B) ∈ Sinc
n of [n], by
distributing n − #B elements of the block [n] \ B among r cycles of σ containing n + 1, . . . , n + r
red elements. Thus,
n
X
#Sn,r = #{(σ1 , B) ∈ Sinc
n : #B = j}r
(n−j)↑

j=0
n   n  
X n (n−j)↑
X n j↑ (n−j)↑
= j!r = 1 r = (1 + r)n↑ ,
j j
j=0 j=0

where the last equality follows from the binomial theorem for raising factorials. By the same
reasoning, for each incomplete permutation (σ1 , B) ∈ Sinc
n ,

#{σ ∈ Sn,r : dperm (σ) = (σ1 , B)} = r(n−#B)↑ .

The proof concludes by taking the ratio. □
3.2. Random r-recursive trees and r-Hoppe trees. The construction is very similar to the
construction of the r-Ewens permutations.
Definition 3.10. Fix n ∈ N0 . An incomplete recursive tree with the root 0 and n additional nodes
is a pair (T, B), where B ⊆ [n] and T is a recursive tree with the root 0 and #B nodes labeled by
the elements of B and colored white. The elements of [n] \ B are called red nodes (these do not
belong to the tree).
Denote by T n the set of all incomplete recursive trees with the root 0 and n additional nodes.
Definition 3.11. Fix parameters r ≥ 0 and τ > 0. A random r-Hoppe tree r-Hoppe[n; τ ] with the
root 0 and n additional nodes is a random element (T, B) of T n with the distribution
r(n−#B)↑ τ degT 0
P[(T, B) = (T, B)] = , (T, B) ∈ T n ,
(r + τ )n↑
where degT 0 is the degree of the root 0 in T . A random r-recursive tree is a special case corre-
sponding to τ = 1.
18
 Note that this definition is correct in a sense that the quantities on the right-hand side sum up
to one. Indeed,
n X j Pn (n−j)↑ n

Pj  j  k
j=0 r k=0 k τ
X r(n−#B)↑ τ deg 0 X r(n−j)↑ τ k X j
n↑
= n↑
1= n↑
(r + τ ) (r + τ ) (r + τ )
(T,B)∈T n j=0 k=0 (T,B)∈T n
#B=j, deg 0=k
Pn n
r(n−j)↑ τ j↑


j=0 j
= = 1,
(r + τ )n↑

where the penultimate equality is a consequence of the fact that there are nj ways to pick j white

nodes and kj ways to construct a recursive tree with the root 0 and j additional white nodes such
 

that the deg 0 = k. The last equality follows from the binomial-type formula for the rising factorial.
The next proposition provides a sequential construction of the r-Hoppe trees.

Proposition 3.12. Consider a Hoppe forest MultiHoppe[n; τ + r; (τ /(τ + r), r/(τ + r))] with d = 2
roots root1 := 0 and root2 having weights θ1 = τ > 0 and θ2 = r ≥ 0, and with n additional nodes.
Let T be the connected component of root1 and B the set of labels of the nodes of T. Then, (T, B)
is distributed as a random tree r-Hoppe[n; τ ].

Proof. Let (T, B) be a fixed incomplete recursive tree. The event {(T, B) = (T, B)} occurs if, and
only if, the nodes with labels in [n] \ B were attached to the subtree rooted at root2 and the nodes
of B formed a particular tree T rooted at root1 . By definition of the Hoppe forest, the probability
of the former event is
Y 1
r(n−#B)↑ ,
j−1+τ +r
j∈[n]\B

whereas the probability of the latter (given the former) is

Y 1
τ degT 0 .
j−1+τ +r
j∈B

Thus, P[(T, B) = (T, B)] = r(n−#B)↑ τ degT 0 /(r + τ )n↑ . The proof is complete. □

Observe that if r = 0, then in the setting of Proposition 3.12 no nodes are attached to root2
(that is, P[#B = n] = 1) and T becomes the usual Hoppe[n; τ ] tree. If, additionally, τ = 1, we
recover RRT[n].
The following is a specialization of the results on multinomial Hoppe forests proved in Section 2.4.

Proposition 3.13. Consider a random r-Hoppe tree (T, B) ∼ r-Hoppe[n; τ ].

τ k n
(a) The root degree of T satisfies deg root1 ∼ r-Stir1[n; τ ], that is, P[deg root1 = k] = (τ +r) n↑ k r
for all k = 0, . . . , n.
(b) The expected number of nodes in T at distance k from the root1 is τ nk τ +r /(τ + r)n↑ .
 

(c) The number of leaves of T, denoted by Tn;1 has the following distribution:
n  
1 X n (n−m)↑
P[Tn;1 = k] = n↑
r A(m − k, k|0, τ ), k = 1, . . . , n.
(τ + r) m
m=k
19
 4. Random r-partitions
Recall that a partition of [n] is a collection B1 , B2 , . . . , Bk of pairwise disjoint, nonempty subsets
Bj ⊆ [n] called blocks, whose union is [n]. The number k of blocks may vary from 1 for the trivial
one-block partition, to n for the finest partition into singletons. Let Πn be the set of all partitions
of [n].
A natural way to generate a partition of [n] is via an urn scheme. Let N ∈ N be a fixed integer.
Drop n balls labeled by the elements of [n] into N urns uniformly at random and independently.
Define a random partition λn,N of [n] by declaring two balls to be in the same block if, and
only if, they fall into the same urn. As has already been mentioned in the introduction, the
number # blocks(λn,N ) of blocks in the resulting partition λn,N has the Stirling–Sibuya distribution
StirSibuya[n; N ] defined by Eq. (1.8):
  k↓
n N
P[# blocks(λn,N ) = k] = = StirSibuya[n; N ]({k}), k ∈ {1, . . . , min(n, N )}. (4.1)
k Nn
4.1. Incomplete partitions derived from urn schemes. Next, we introduce an r-version
of (4.1). To this end, it is necessary to replace partitions by incomplete partitions.

Definition 4.1. An incomplete partition of [n] is a collection (B0 , {B1 , B2 , . . . , Bk }) of pairwise

disjoint subsets Bj ⊆ [n], whose union is [n] and Bj ̸= ∅ for all j ∈ {1, . . . , k}. Note that B0 , called
the red block, is allowed to be empty. The blocks B1 , . . . , Bk , called the white blocks, are non-empty.
Let Πinc
n be the set of all incomplete partitions of [n].

The already constructed r-versions of Ewens permutations and random recursive trees suggest
the following way to generate a random incomplete partition leading to an r-StirSibuya[n; N ]-
distribution. We consider an urn scheme with N + 1 urns 0, 1, . . . , N . The urn 0 is red and has
frequency r/(N + r) ≥ 0, whereas urns 1, 2, . . . , N are white and have frequencies 1/(N + r). Drop
n balls into the urns and define an incomplete partition λn,N,r of [n] by declaring the balls in the
urn 0 to form the red block B0 whereas the balls in those urns 1, . . . , N that are non-empty to form
white blocks as before.

Proposition 4.2. Let N ∈ N be an integer and r ≥ 0. Then,

r#B0 N k↓
P[λn,N,r = (B0 , {B1 , . . . , Bk })] = , (B0 , {B1 , . . . , Bk }) ∈ Πinc
n . (4.2)
(N + r)n
Proof. The formula follows upon noticing that
#B0 
r#B0 N k↓ N k↓
 
r
= .
(N + r)n N +r (N + r)n−#B0
The first factor on the right-hand side is the probability for balls with labels in B0 to fall in the
urn 0. The second factor is the probability for remaining balls to form a partition B1 , . . . , Bk of
[n] \ B0 . □

Remark 4.3. For r = 0, the right-hand side of (4.2) vanishes whenever #B0 ̸= 0. Thus, neglecting
the (empty) red block in λn,N,0 we obtain the usual partition λn,N of [n].

Let type be the mapping which sends an incomplete partition π ∈ Πinc

n to the multiset type(π)
of block sizes of π with a distinguished part b0 ∈ N0 . Formally, an incomplete partition π =
20
(B0 , {B1 , . . . , Bk }) has type type(π) = (b0 , {b1 , . . . , bk }), where b0 ∈ N0 and b1 , . . . , bk ∈ N if, and
only if, #B0 = b0 and {b1 , . . . , bk } is the multiset of block sizes of B1 , . . . , Bk .
From Proposition 4.2 we immediately conclude the following.

Corollary 4.4. Let N ∈ N be an integer and r ≥ 0. Then,

1 rb0 N k↓
 
n
P[type(λn,N,r ) = (b0 , {b1 , . . . , bk })] = (4.3)
b0 , . . . , bk k! (N + r)n
for each collection (b0 , {b1 , . . . , bk }) of integers such that b0 ∈ N0 , b1 , . . . , bk ∈ N and b0 + b1 + · · · +
bk = n.

For an incomplete partition π ∈ Πinc

n let # blocks(π) denote the number of white blocks.

Proposition 4.5. Let N, n ∈ N be integers and r ≥ 0. The number of white blocks in λn,N,r has
the r-StirSibuya[n; N ] distribution, that is,
N k↓
 
n
P[# blocks(λn,N,r ) = k] = , k ∈ {0, 1, . . . , min(n, N )}. (4.4)
k r (N + r)n
Similarly to Proposition 3.5 we give three simple proofs of this assertion exploiting various
representations of nk r .


First proof. This proof uses the second equality in (3.4). One just needs to sum the right-hand sides
of (4.3) over all (k +1)-tuples (b0 , b1 , . . . , bk ) such that b0 ∈ N0 , b1 , . . . , bk ∈ N and b0 +b1 +· · ·+bk =
n. Formula (4.4) follows then from (3.4). □
Pn
Second proof. This proof hinges on an equality in (1.6), namely, k r = j=k nj kj rn−j . Condi-
n


tioning on the number of balls that do not fall in the urn 0 and using (4.1), we obtain

P[# blocks(λn,N,r ) = k] = StirSibuya[Bin[n; N/(N + r)]; N ]({k})

n   n−j  j   k↓
N k↓
 
X n r N j N (1.6) n
= =
j N +r N +r k Nj k r (N + r)n
j=k

for each k ∈ {0, 1, . . . , min(n, N )}. □

n
Third proof. In this proof we suppose that r ∈ N is integer and exploit the formula k r =
Pn n j
(j−k)↓
j=k j k r , see (1.6). Replace the red urn with frequency r/(N + r) by r urns with
frequency 1/(N + r), so that we allocate n balls in N + r equiprobable urns (N white ones and r
red ones). In total, there are (N + r)n possible allocations. Then,
n
X
P[# blocks(λn,N,r ) = k] = P[k white urns and j − k red urns are non-empty]
j=k
n   
N k↓
 
1 X n j k↓ (j−k)↓ n
= N r = ,
(N + r)n j k k r (N + r)n
j=k

where j is the number of ways to decompose n balls into j blocks, kj is the number of ways to
n

choose k blocks to be placed into white urns, and, finally, N k↓ , respectively r(j−k)↓ , is the number
of ways to choose white, respectively red, urns to be filled. □

Remark 4.6. Proposition 4.5 can also be found in [53, Section 5].
21
4.2. Random uniform r-partitions and random Gibbs r-partitions. The cardinality of the
set Πn of all partitions of [n] is the Bell number Bn . Let π be a random uniform partition of [n],
that is, a random element with values in Πn and such that P[π = π] = 1/Bn for each partition
π ∈ Πn . Stam [63] discovered a randomization of the number of urns N in the above scheme which
leads to a random partition with the uniform distribution on Πn . Let M be a random variable with
1 mn
P[M = m] = , m ∈ N.
e · Bn m!
The numbers on the right-hand side sum up to 1 by the Dobiński formula. Stam showed that the
distribution of λn,M is uniform on Πn , where M is assumed independent of everything else. The
proof is rather simple. From (4.2) applied with r = 0 we conclude that, for each {B1 , . . . , Bk } ∈ Πn ,
∞ ∞
1 X mn mk↓ 1 X 1 1
P[λn,M = {B1 , . . . , Bk }] = = = .
e · Bn m! mn e · Bn (m − k)! Bn
m=k m=k

Moreover, Stam showed that the number of empty boxes is Poisson distributed with parameter 1
and independent of the partition.
We now focus on a two-fold generalization of Stam’s construction. This leads us to the r-Stirling
distributions of the second kind r-Stir2[n; θ].

Definition 4.7. Fix some θ > 0 and r ≥ 0. A random variable πθ,r with values in Πinc
n is called a
random Gibbs r-partition if it has the distribution
θk r#B0
P[πθ,r = (B0 , {B1 , . . . , Bk })] = , (B0 , {B1 , . . . , Bk }) ∈ Πinc
n ,
Tn,r (θ)
where we recall the definition
n  
X n
Tn,r (θ) = θk
k r
k=0
of the r-Touchard polynomial. This distribution is denoted by r-GibbsPart[n; θ]. If θ = 1, an
incomplete partition π1,r is called uniform r-partition. The number Tn,r (1) =: Bn,r is called the
n-th r-Bell number, see [50, § 8.3] and [48].

Remark 4.8. In the literature, see, for example, [57] and [19], it is more common to use the term
‘Gibbs partition’ in a more general sense when the weight of a block also depends on its size. In
Definition 4.7 we use this terminology in a narrow sense, when all white blocks have the same (but
arbitrary) weight θ > 0.

Note that the above distribution resembles (4.2) with N replaced by θ. The change of the factor
θk↓ to θk results in a more sophisticated normalizing factor Tn,r (θ) in place of (θ + r) n in (4.2).

The correctness of the definition follows from the chain of equalities (there are nj choices for a

red block B0 containing n − j elements and nk choices to partition the remaining j elements into


k white blocks B1 , . . . , Bk ):
j
n X    n n   
X 1 X
k n−j n j 1 X
k
X
n−j n j
P[πθ,r = π] = θ r = θ r
Tn,r (θ) j k Tn,r (θ) j k
π∈Πninc j=0 k=0 k=0 j=k
n  
(1.6) 1 X n
= θk = 1.
Tn,r (θ) k r
k=0
22
Proposition 4.9. Let πθ,r ∼ r-GibbsPart[n; θ]. Then

rb0 θk n! rb0 θk
 
n 1 1
P[type(πθ,r ) = (b0 , {b1 , . . . , bk })] = = . (4.5)
Tn,r (θ) b0 , . . . , bk k! k! b0 ! · · · bk ! Tn,r (θ)

Furthermore, the number of white blocks follows the r-Stirling distribution of the second kind:

θk
 
n
P[# blocks(πθ,r ) = k] = , k = 0, . . . , n. (4.6)
Tn,r (θ) k r

Proof. Equality (4.5) is obvious. Formula (4.6) follows from the second equality in (3.4) upon
summing the right-hand side of (4.3) over all b0 ∈ N0 and b1 , . . . , bk ∈ N satisfying b0 + · · · + bk =
n. □

Next, in the spirit of Stam’s construction, we show that πθ,r can be realized by randomizing
N in the incomplete partition λn,N,r defined at the beginning of this section. Consider a random
variable Mθ,r with the following distribution:

1 (r + m)n θm
P[Mθ,r = m] = , m ∈ N0 . (4.7)
eθ Tn,r (θ) m!

We first check that this is indeed a probability distribution. The next lemma is an r-analogue
of the Dobiński formula, see [48, Theorem 5.1]. For completeness, we provide a proof.

Lemma 4.10. For each θ > 0 and r ≥ 0,

∞
X (r + m)n θm
= eθ Tn,r (θ).
m!
m=0

Proof. Using the second formula in (1.2) yields

∞ ∞ n   n   X ∞
X (r + m)n θm X θm X n X n θm k↓
= mk↓ = m
m! m! k r k r m!
m=0 m=0 k=0 k=0 m=0
n   ∞ n   ∞
! !
X n X θ m−k X n X θm
= θk = θk = Tn,r (θ)eθ .
k r (m − k)! k r m!
k=0 m=k k=0 m=0

Here is a version of Stam’s result including the claim that the number of empty urns is in-
dependent of the incomplete partition and has a Poisson distribution. Naturally, Mθ,r is taken
independent of everything else.

Proposition 4.11. For each θ > 0 and r ≥ 0,

λn,Mθ,r ,r , Mθ,r − # blocks(λn,Mθ,r ,r ) ∼ r-GibbsPart[n; θ] ⊗ Poi[θ].
23
Proof. Fix v ∈ N0 and an incomplete partition (B0 , {B1 . . . , Bk }) ∈ Πinc
n . Then, by (4.7) and
Proposition 4.2,
P[λn,Mθ,r ,r = (B0 , {B1 . . . , Bk }), Mθ,r = k + v]
= P[Mθ,r = k + v] · P[λn,Mθ,r ,r = (B0 , {B1 . . . , Bk })|Mθ,r = k + v]
= P[Mθ,r = k + v] · P[λn,k+v,r = (B0 , {B1 . . . , Bk })]
1 (r + k + v)n θk+v r#B0 (k + v)k↓
= ·
eθ Tn,r (θ) (k + v)! (k + v + r)n
θk r#B0 −θ θv
= ·e
Tn,r (θ) v!
= r-GibbsPart[n; θ]((B0 , {B1 . . . , Bk })) · Poi[θ]({v}),
where in the last line we have used Definition 4.7. □
As we did for r-permutations, we conclude the discussion of r-partitions with a remark on
Broder’s definition of the r-Stirling numbers of the second kind given in the introduction. If r ∈ N0
is an integer, then the r-uniform partition can be constructed combinatorially as follows. Consider
the set [n + r] in which elements 1, 2, . . . , n are white and n + 1, . . . , n + r are red. Let Πn,r be the set
of all partitions of [n + r] such that red elements are in different blocks and let dpart : Πn,r 7→ Πinc n
be the mapping defined as follows. For π ∈ Πn,r , the image dpart (π) is an incomplete partition
(Bb0 , {B bk }) of [n] such that the red block B
b1 , . . . , B b0 is obtained by merging r blocks of π containing
red elements, removing these red elements from the union, and keeping white blocks B b1 , . . . , B
bk
unchanged. Let π bn,r be a partition of [n + r] picked uniformly at random from Πn,r .

Proposition 4.12. For each n ∈ N and r ∈ N0 ,

πn,r ) ∼ r-GibbsPart[n; 1],
dpart (b
that is, dpart (b
πn,r ) is the uniform r-partition.

Proof. Each π ∈ Πn,r can be constructed from an incomplete partition π1 := (B0 , {B1 , . . . , Bk }) ∈
Πinc
n of [n] by distributing b0 = #B0 elements of the block B0 among r blocks of π containing r red
elements (there is an r-fold choice, for each x ∈ B0 ). Thus,
n n  
X
inc j
X n
#Πn,r = #{π1 ∈ Πn : #B0 = j}r = Bn−j rj = Bn,r = Tn,r (1),
j
j=0 j=0

where the penultimate inequality follows from (1.6). Furthermore, for each incomplete partition
π1 ∈ Πinc
n with red block B0 of size b0 ,

#{π ∈ Πn,r : dpart (π) = π1 } = rb0 .

The proof concludes by taking the ratio. □

5. Random r-compositions
A composition of n into k summands is a tuple (b1 , . . . , bk ) with b1 , . . . , bk ∈ N and b1 + · · · + bk =
n. Using stars-and-bars argument, one can check that the number of compositions of n into k
n−1


summands is k−1 . Picking one of these compositions at random, we obtain a uniform random
composition of n into k summands. In the following we introduce an r-analogue of this notion.
24
Definition 5.1. An incomplete composition of n is a tuple (b0 , b1 , . . . , bk ) such that b0 ∈ N0 ,
b1 , . . . , bk ∈ N and b0 + b1 + · · · + bk = n. Note that b0 is allowed to be 0, while b1 , . . . , bk are not.
The number of such incomplete compositions is nk by stars-and-bars argument. Incomplete

compositions are closely related to Weyl chambers of type B. Consider Weyl chambers of types A
and B which are polyhedral cones in Rn given by A(n) := {(x1 , . . . , xn ) ∈ Rn : x1 ≤ · · · ≤ xn } and
B (n) := {(x1 , . . . , xn ) ∈ Rn : 0 ≤ x1 ≤ · · · ≤ xn }. Then, the k-dimensional faces of the cone A(n)
are in bijective correspondence with compositions of n into k summands and have the form
{x1 = · · · = xb1 ≤ xb1 +1 = · · · = xb1 +b2 ≤ · · · ≤ xb1 +···+bk−1 +1 = · · · = xn }.
Similarly, k-dimensional faces of B (n) are in bijective correspondence with the incomplete compo-
sitions (b0 , b1 , . . . , bk ) and have the form

{0 = x1 = · · · = xb0 ≤ xb0 +1 = · · · = xb0 +b1 ≤ xb0 +b1 +1 = . . .

= xb0 +b1 +b2 ≤ · · · ≤ xb0 +b1 +···+bk−1 +1 = · · · = xn }.
(n,k) (n,k) (n,k)
Definition 5.2. Fix r ≥ 0. A random r-composition of n is a random vector (b0 , b1 , . . . , bk )
such that, for each incomplete composition (b0 , b1 , . . . , bk ) of n,
b0 +r−1


(n,k) (n,k) (n,k) rb0 ↑ /b0 ! b0
P[b0 = b0 , b1 = b1 , . . . , bk = bk ] = n+r−1
= n+r−1

. (5.1)
k+r−1 k+r−1
n−1


Example 5.3. For r = 0, we understand the right-hand side as 1/ k−1 for b0 = 0 and as 0 for
(n,k)
b0 ≥ 1. So, b0 = 0 with probability 1 and, discarding the 0-th entry, we recover the uniform
distribution on the set of compositions of n.
Example 5.4. For r = 1, each incomplete composition (b0 , b1 , . . . , bk ) is equally likely, with the
n


corresponding probability being 1/ k . Thus, we recover the uniform distribution on the set of
all incomplete compositions which are in bijective correspondence with the k-faces of the Weyl
chamber B (n) . Note in passing that it would be possible to give a similar interpretation for the
r-uniform partitions with r = 1/2 by considering uniformly distributed elements of the subspace
lattice generated by the type B reflection arrangement, see [1, § 6.3, 6.7] for a description of the
correspondence between (incomplete) partitions and flats in reflection arrangements.
With the help of the next lemma we will show that Definition 5.2 indeed defines a probability
distribution.
Lemma 5.5. Let r ≥ 0. Then, for each m ∈ N and ℓ ∈ {0, 1, . . . , m},
m−ℓ
X b0 + r − 1m − b0  m−ℓ X b0 + r − 1m − b0  m + r
= = . (5.2)
b0 ℓ r−1 ℓ ℓ+r
b0 =0 b0 =0

Proof. The first equality is trivial and the second one is Eq. (5.26) in [22], provided that r ∈ N.
The latter constraint is superfluous (indeed, both sides are polynomials in r). For the reader’s
convenience we give below yet another proof of (5.2) valid for all real r ≥ 0. We use the Taylor
expansions
∞  ∞  
xℓ

X b0 + r − 1 b0 −r
X j j
x = (1 − x) and x = .
b0 ℓ (1 − x)ℓ+1
b0 =0 j=ℓ
25
Multiplying these Taylor series and evaluating the coefficient of xm , we arrive at the left-hand side
of (5.2). On the other hand, the same coefficient can be computed as follows:
∞
(1 − x)−r xℓ xj (ℓ + r + 1)j↑
   
m m−ℓ −ℓ−r−1 m−ℓ
X
j↑ m+r
[x ] = [x ](1−x) = [x ] (ℓ+r+1) = = .
(1 − x)ℓ+1 j! (m − ℓ)! ℓ+r
j=0

Comparing these results completes the proof. □

Remark 5.6. Now we can show that (5.1) indeed defines a probability distribution on the set
of incomplete compositions of
n. The number of incomplete compositions with a given value of
0 −1
b0 ∈ {0, . . . , n − k} is n−b
k−1 . Each such an incomplete composition has probability given by the
right-hand side of (5.1). The sum of all such probabilities is
n−k
X n − b0 − 1b0 + r − 1
1
n+r−1 =1
k−1


k+r−1
b0
b0 =0
by Lemma 5.5 with m = n − 1 and ℓ = k − 1.
5.1. Couplings of random r-compositions and the Dirichlet multinomial distribution. In
this subsection we aim at constructing an r-composition via a certain urn scheme, in a way similar
to that used for uniform r-partitions in Section 4. As a result, we provide several couplings of
r-compositions which are consistent either in n or in k. It turns out that the random r-composition
can be constructed using an urn scheme with random frequencies.
(n,k) (n,k) (n,k)
Formula (5.1) implies that the random vector (b0 , b1 − 1, . . . , bk − 1) has the Dirichlet
multinomial distribution
DirMult[n − k; r, 1, . . . , 1], (5.3)
| {z }
k times
see definition (2.7). It is known that the probability measure DirMult[n; α0 , α1 , . . . , αk ] can be
regarded as a mixture of multinomial distributions. More precisely,
DirMult[n; α0 , α1 , . . . , αk ] = Mult[n; Dir[α0 , α1 , . . . , αk ]],
where the Dirichlet distribution Dir[α0 , α1 , . . . , αk ] is defined by the density
k
Γ(α0 )Γ(α1 ) · · · Γ(αk ) α0 −1 α1 −1
· · · xαk k −1 ,
X
(x0 , . . . , xk ) 7→ x x1 xi = 1, xi > 0. (5.4)
Γ(α0 + α1 + · · · + αk ) 0
i=0
The above interpretation suggests a coupling of random r-compositions which is consistent in
n ≥ k if k is fixed. The coupling works as follows. Take a random partition P of [0, 1) into k + 1
intervals (urns)
Ik,0 := [0, P0 ), Ik,1 := [P0 , P0 + P1 ), ... Ik,k := [P0 + · · · + Pk−1 , 1),
where (P0 , P1 , . . . , Pk ) follows Dir[r, 1, . . . , 1] distribution. In a role of balls take an independent of
P sample U1 , U2 , . . . , Un , . . . from the uniform distribution on [0, 1]. Let Zn,i , i = 0, . . . , k be the
number of balls among U1 , . . . , Un in the urn Ik,i , that is,
Zn,i := #{1 ≤ j ≤ n : Uj ∈ Ik,i }, 0 ≤ i ≤ k.
The following trivially holds true.
Proposition 5.7. Fix k ∈ N. For each n ≥ k, the random vectors (Zn−k,0 , Zn−k,1 , . . . , Zn−k,k ) and
(n,k) (n,k) (n,k)
(b0 , b1 − 1, . . . , bk − 1) have the same Dirichlet multinomial distribution (5.3).
26
 The aggregation property of the Dirichlet distribution asserts that if
(X0 , X1 , . . . , Xk ) ∼ Dir[α0 , α1 , . . . , αk ],
then
(X0 + X1 , X2 , . . . , Xk ) ∼ Dir[α0 + α1 , α2 , . . . , αk ].
This property suggests another coupling of r-compositions which is consistent in both parameters
n and k but works only when r is an integer. This coupling hinges on a simple observation that the
collection of gaps between consecutive order statistics from the uniform distribution on [0, 1] has
a symmetric Dirichlet distribution Dir[1, 1, . . . , 1]. More precisely, let V1 , V2 , . . . be a sample from
the uniform distribution on (0, 1). For each M ∈ N, let
0 < VM,1 < VM,2 < · · · < VM,M < 1
be the order statistics of the first M points of the sample. Then
(VM,1 , VM,2 − VM,1 , . . . , VM,M − VM,M −1 , 1 − VM,M ) ∼ Dir[ 1, . . . , 1 ].
| {z }
M +1 times
By the aggregation property, if r is an integer
(Vk+r,r , Vk+r,r+1 − Vk+r,r , . . . , Vk+r,k+r − Vk+r,k+r−1 , 1 − Vk+r,k+r ) ∼ Dir[r, 1, . . . , 1 ].
| {z }
k+1 times

Thus, a family (Ik,0 , Ik,1 , . . . , Ik,k ), k ∈ N of nested partitions of [0, 1) such that the lengths of
(Ik,0 , Ik,1 , . . . , Ik,k ) have the Dir[r, 1, . . . , 1] distribution can be sequentially constructed as follows.
Start with a sample V1 , . . . , Vr of size r from the uniform distribution on [0, 1]. This sample yields
a partition of [0, 1) into two intervals I1,0 := [0, Vr,r ) and I1,1 := [Vr,r , 1), whose lengths follow the
Dir[r, 1] distribution. Add a new point Vr+1 . If Vr+1 lands into I1,1 , then it splits I1,1 into two
subintervals I2,1 := [Vr,r , Vr+1 ) and I2,2 := [Vr+1 , 1). In this case, put also I2,0 := [0, Vr,r ). On
the other hand, if Vr+1 lands into I1,0 , the splitting is I2,0 := [0, Vr+1,r ) and I2,1 := [Vr+1,r , Vr,r ).
In this case, put also I2,2 := I1,1 . More generally, if the partition (Ik,0 , Ik,1 , . . . , Ik,k ) has already
been constructed, adding a point Vr+k yields a partition (Ik+1,0 , Ik+1,1 , . . . , Ik+1,k+1 ) depending on
whether Vr+k ∈ Ik,0 or Vr+k ∈ / Ik,0 . Formally, the construction is given by the equations

Ik+1,0 := [0, Vr+k,r ), Ik+1,j := [Vr+k,r+j−1 , Vr+k,r+j ), j = 1, . . . , k,

Ik+1,k+1 = [Vr+k,r+k , 1). (5.5)
We used the word ‘nested’ to highlight that the partition (Ik+1,0 , Ik+1,1 , . . . , Ik+1,k+1 ) is obtained
from the partition (Ik,0 , Ik,1 , . . . , Ik,k ) by splitting one of the intervals Ik,j into two subintervals. By
taking samples (Vi ) and (Ui ) independent and appealing to Proposition 5.7 we obtain (in the case
r ∈ N) a coupling of r-compositions with the following properties:
(n,k+1) (n,k+1) (n,k+1) (n,k) (n,k) (n,k)
• an r-partition (b0 , b1 , . . . , bk+1 ) is derived from an r-partition (b0 , b1 , . . . , bk )
(n,k)
by splitting one of blocks bj into two parts;
(n+1,k) (n+1,k) (n+1,k) (n,k) (n,k) (n,k)
• an r-partition (b0 , b1 , . . . , bk ) is derived from (b0 , b1 , . . . , bk ) by increasing
one of blocks by one without changing other blocks.
(n,k) (n,k) (n,k)
Proposition 5.7 is a de Finetti-type result, which provides a realization of (b0 , b1 , . . . , bk )
via a randomized multinomial allocation scheme. Given next is a Markovian construction that
(n,k) (n,k) (n,k)
gives a realization of the process (b0 , b1 , . . . , bk )n≥k , with k fixed, as consecutive values of
a certain Markov chain.
27
Proposition 5.8. Fix k ∈ N. Let (X(l))l∈N0 be a Markov chain on [0, ∞) × Nk with the initial
state X(0) = (r, 1, . . . , 1) and transition probabilities
ij
P[X(l + 1) = (i0 , i1 , . . . , ik ) + ej |X(l) = (i0 , i1 , . . . , ik )] = ,
i0 + i1 + · · · + ik
j = 0, . . . , k, l ∈ N0 ,
where ej is the j-th unit vector in Rk+1 . For each fixed n ≥ k, X(n − k) − re0 has the same
(n,k) (n,k) (n,k)
distribution as (b0 , b1 , . . . , bk ).
The Markov chain constructed in Proposition 5.8 is nothing else but a Pólya’s urn. Initially, the
urn contains r balls of color 0 and one ball of each of other colors 1, 2, . . . , k. At each step, a ball
is drawn uniformly at random from the urn, its color is observed, and the ball is returned to the
urn with an additional ball of the same color. The vector X(n) encodes the colors of balls in the
urn after n steps.
For later needs, we formulate a lemma.
(n) (n)
Lemma 5.9. Fix n, k ∈ N. The random vector (X0 , . . . , Xk ) with the DirMult[n; α0 , α1 , . . . , αk ]
distribution is negatively associated, that is, for each pair of disjoint subsets A1 , A2 of {0, 1, . . . , k}
and all nondecreasing in each coordinate functions f : R|A1 | → R and g : R|A2 | → R
(n) (n) (n) (n)
E[f (Xi , i ∈ A1 )g(Xj , j ∈ A2 )] ≤ E[f (Xi , i ∈ A1 )]E[g(Xj , j ∈ A2 )].
Proof. This claim is given on p. 292 in [33] without proof. We now give a proof. Let Y0 , Y1 , . . . , Yk be
independent random variables such that, for i ∈ {0, . . . , k}, Yi has a negative binomial distribution
with parameters αi and β ∈ (0, 1), that is,
αik↑ αi
P[Yi = k] = NBin[αi ; β]({k}) = β (1 − β)k , k ∈ N0 .
k!
(n) (n)
We claim that the distribution of (X0 , . . . , Xk ) coincides with the conditional distribution of
(Y0 , . . . , Yk ) given ki=0 Yi = n. Indeed, for each (b0 , b1 , . . . , bk ) ∈ Nk+1
P
0 summing up to n,
k
" #
X P[Y0 = b0 ]P[Y1 = b1 ] · · · P[Yk = bk ]
P Y0 = b0 , Y1 = b1 , . . . , Yk = bk Yi = n =
P[ ki=0 Yi = n]
P
i=0
k
!
Y
= NBin[αi ; β]({bi }) /NBin[α0 + · · · + αk ; β]({n}).
i=0
After cancellations this becomes the right-hand side of (2.7) with d = k + 1, ki = bi and θi = αi ,
for i = 0, . . . , k. In view of the above representation we may use [33, Theorem 2.6], which states it
is sufficient to show that for each subset A of {0, 1, . . . , k} and all nondecreasing in each coordinate
functions f : R|A| → R the sequence
h k i
(n)
X
n 7→ E f (Yi , i ∈ A) Yj = n = E[f (Xi , i ∈ A)]
j=0

is nondecreasing. This property is secured by the aforementioned coupling, which enables us to

(n)
think of Xi as the occupancy count of the box i corresponding to (P0 , P1 , . . . , Pk ) having the
Dir[α0 , α1 , . . . , αk ] distribution. Throwing a new ball does not decrease an occupancy count, that
(n) (n+1)
is, Xi ≤ Xi a.s. □
28
 Lemma 5.9 in combination with Proposition 5.7 immediately implies a result to be used later in
the proof of Lemma 6.23.
(n,k) (n,k) (n,k)
Corollary 5.10. The random r-composition (b0 , b1 , . . . , bk ) is negatively associated.
5.2. Marginal distributions of r-compositions and their limits.
(n,k) (n,k) (n,k)
Proposition 5.11. The marginal distributions of the random r-composition (b0 , b1 , . . . , bk )
are given by
n−b0 −1 b0 +r−1



h i
(n,k) k−1 b0
P b0 = b0 = n+r−1

, b0 ∈ {0, 1, . . . , n − k}, (5.6)
k+r−1
h i n−bj +r−1

(n,k) k+r−2
P bj = bj = n+r−1

, bj ∈ {1, . . . , n − k + 1}. (5.7)
k+r−1
Moreover, for 0 < i < j ≤ k, the bivariate distributions are given by
h i n+r−bi −bj −1

(n,k) (n,k) k+r−3
P bi = bi , bj = bj = n+r−1

, bi , bj ∈ {1, . . . , n − k + 1}, bi + bj ≤ n − k + 2.
k+r−1
(5.8)
n−b0 −1


Proof. For (5.6), just observe that the number of incomplete compositions with fixed b0 is k−1
and the probability of each such an incomplete composition is given by the right-hand side of (5.1).
For (5.7), we can take j = 1 without loss of generality. In the case k = 1, (5.6) secures
n−b1 +r−1 n−b1 +r−1



h i h i
(n,k) (n,k) n−b1 r−1
P b1 = b1 = P b0 = n − b1 = n+r−1

= n+r−1

.
k+r−1 k+r−1
Assume now that k ≥ 2.
Observe that the number of incomplete compositions with fixed b0 and
0 −b1 −1
b1 is given by n−bk−2 . Since the probability of each such an incomplete composition is given
by the right-hand side of (5.1), we obtain
i n−bX 1 −k+1   b0 +r−1
n−b1 −1+r


h
(n,k) n − b0 − b1 − 1 b0 k−2+r
P b1 = b1 = n+r−1 = n+r−1

,
k−2


b0 =0 k+r−1 k+r−1

where we evaluated the sum using Lemma 5.5 with m := n − b1 − 1 and ℓ := k − 2.

Alternatively, one can use a known property of the Dirichlet distribution Dir[α0 , α1 , . . . , αk ],
which asserts that its i-th component has the beta distribution Beta[αi ; α−αi ], where α := kj=0 αj .
P

Then, by Proposition 5.7,

(n,k) (n,k)
b0 ∼ Bin[n − k; Beta[r; k]] and bj ∼ 1 + Bin[n − k; Beta[1; k + r − 1]]. (5.9)
Formulas (5.6) and (5.7) follow by simple manipulations with beta-integrals.
For the proof of (5.8), note that by the aggregation property of the Dirichlet distribution in
combination with Proposition 5.7,
 
b0(n,k) + (n,k) (n,k) (n,k)
X
(bℓ − 1), bi − 1, bj − 1 ∼ DirMult[n − k; Dir[r + k − 2, 1, 1]]
ℓ∈{0,i,j}
/

for each pair of indices 0 < i < j ≤ k. Thus,

ZZZ
(n,k)
bi
(n,k)
−1 bj −1 Γ(r + k)
E[t s ]= (x + ty + sz)n−k xr+k−3 dxdydz, t, s ∈ C,
Γ(r + k − 2)
29
where the integration is taken over the simplex x + y + z = 1, x, y, z ≥ 0. By expanding (x + ty +
sz)n−k with the aid of the multinomial theorem and equating the coefficients, we conclude that
h
(n,k) (n,k)
i Γ(r + k) (n − k)! Γ(n + r − bi − bj )
P bi = bi , bj = bj = .
Γ(r + k − 2) (n − k + 2 − bi − bj )! Γ(n + r)
This yields (5.8) after elementary manipulations. □
(n,k) (n,k) (n,k)
Corollary 5.12. The conditional distribution of (b1 , . . . , bk ) given b0 is uniform on the
(n,k)
set of all (usual) compositions of n − b0 into k blocks.
(n,k) j −1
= bj ] = n−b

n−1

Remark 5.13. For r = 0, P[bj k−2
/ k−1 for all j ∈ {1, . . . , k}. Now observe that
(n,k)
the formula for the marginal distribution of bj for general r ≥ 0 can be obtained from the r = 0
case if we replace n and k by n + r and k + r, respectively. This observation has several implications
given below.
Proposition 5.14. For r ≥ 0,
(n,k) n+r (n,k) r(n − k)
E[bj ]= , j ∈ {1, . . . , k}, E[b0 ]= . (5.10)
k+r k+r
(n,k)
Proof. For r = 0 it is clear, by exchangeability, that E[bj ] = n/k for all j ∈ {1, . . . , k}. According
to Remark 5.13, for general r ≥ 0, we have to replace n and k by n + r and k + r, which gives the
(n,k) (n,k) (n,k)
first formula in (5.10). The second formula follows from b0 = n − b1 − · · · − bk . □
A strong law of large numbers for the multinomial distribution immediately implies a limit
theorem for r-compositions with k being fixed.
Proposition 5.15. Let k ∈ N be a fixed integer and r ≥ 0. Then
1 (n,k) (n,k) (n,k) d
(b0 , b1 , . . . , bk ) −→ Dir[r, 1, . . . , 1].
n n→∞ | {z }
k times

Proposition 5.16. Fix some r ≥ 0 and let n → ∞, k → ∞ such that k/n → α for some α ∈ (0, 1].
(n,k) (n,k) (n,k)
Then, for a random r-composition (b0 , b1 , . . . , bk ),
h i b + r − 1
(n,k) 0
lim P b0 = b0 = (1 − α)b0 αr , b0 ∈ {0, 1, . . . },
n→∞ b0
h i
(n,k)
lim P bj = b1 = α(1 − α)b1 −1 , b1 ∈ {1, 2, . . . }, j ∈ {1, . . . , k}.
n→∞
(n,k) (n,k)
That is to say, b0 converges weakly to the negative binomial distribution NBin[r; α], while bj
converges weakly to the geometric distribution Geom[α].
Proof. Follows from Proposition 5.11 after letting n → ∞. □
Proposition 5.17. Fix some r ≥ 0 and let n → ∞, k → ∞ such that k/n → 0. Then, for a
(n,k) (n,k) (n,k)
random r-composition (b0 , b1 , . . . , bk ),
(n,k) (n,k)
b0 d b1 d
−→ Gamma[r; 1] and −→ Exp[1], (5.11)
n/k n,k→∞ n/k n,k→∞

where Gamma[r; 1] is the gamma distribution with parameters r and 1 and Exp[1] is the standard
exponential distribution.
30
Proof. The easiest way to prove this is to use representations (5.9). The first claim in (5.11) follows
from the next three facts. First, with (γ(t))t≥0 being the standard gamma-subordinator,
γ(r)
∼ Beta[r; k],
γ(r + k)
see [41, Section 9.1]. Second, for each fixed λ > 0,
k P
Bin[n − k; λ/k] → λ, k, n → ∞, k/n → 0,
n
as can be checked by using Chebyshev’s inequality. Third, by the strong law of large numbers
γ(r + k) a.s.
−→ 1.
k k→∞

It remains to note that γ(r) ∼ Gamma[r; 1].

The second claim in (5.11) follows by the same reasoning with the help of the second relation
in (5.9). The proof is complete. □

6. The (r, s)-Lah distribution

6.1. Motivation. In this section we investigate the so-called (r, s)-Lah distribution which is defined
in terms of the r-Stirling numbers of both kinds. Special cases of this distribution appear in [24] in
the context of random recursive trees (and the closely related shortest path trees on the complete
graph with exponentially distributed weights), and in [36, 38]. In the last two articles, limit
theorems for this distribution are applied to obtain asymptotics of the face numbers of random walk
convex (and positive) hulls [20, 21, 39, 40]. All these papers investigate the (r, s)-Lah distribution
for certain special values of parameters, for example, the distribution appearing in [24] corresponds
to (r, s) = (1, 0), the distributions in [36, 39, 40] and [21] correspond to (0, 0) and (1/2, 1/2),
respectively, while [38] is concerned with the case r = s.
Our purpose is to introduce a general distributional family unifying all these special cases. We
relate these general (r, s)-Lah distributions to the (r + s)-versions of the classic combinatorial
probability structures defined in the previous sections. In particular, we derive several stochastic
representations of the (r, s)-Lah distribution in terms of random (r + s)-compositions and (multi-
nomial) Hoppe trees. In Section 7 we use these representations when proving limit theorems for
the (r, s)-Lah distribution.

6.2. The r-Lah numbers. We start by reviewing some properties of the r-Lah numbers L(n, k)r ;
see [54] as well as [3, 14, 59] for more information. For n ∈ N0 , k ∈ {0, . . . , n} and r ∈ R, these
numbers may be defined [54, Theorem 3.10] by their exponential generating function
∞ k
xn

X 1 x
L(n, k)r = (1 − x)−2r . (6.1)
n! k! 1 − x
n=k

We extend this definition by putting L(n, k)r := 0 for n ∈ N0 and k ∈

/ {0, . . . , n}. The Taylor
expansion of the right-hand side of (6.1) gives [54, Theorem 3.7]
 
n + 2r − 1 n!
L(n, k)r = . (6.2)
k + 2r − 1 k!
31
As a consequence, L(n, k)r is a polynomial of r, for fixed n and k. The following formula [54,
Theorem 3.2] complements (1.2):
n
X
(x + r)n↑ = L(n, k)r (x − r)k↓ . (6.3)
k=0

Using (1.2) and (6.3) it can be shown [54, Theorem 3.11] that
n    
X n j
L(n, k)r = . (6.4)
j r k r
j=k

More generally, from the same reference it is known that

n    
X n j
L(n, k) r+s = . (6.5)
2 j r k s
j=k

The cited paper [54] only considers nonnegative integer r and requires that r and s have the same
parity. In fact, since both sides of the identity are polynomials in r and s, it holds for arbitrary
r, s ∈ R.
From Eq. (5.1) in the definition of r-compositions we conclude that
n! X (2r)ℓ0 ↑
L(n, k)r = ,
k! ℓ0 !
ℓ0 ∈N0 ,ℓ1 ,...,ℓk ∈N
ℓ0 +ℓ1 +···+ℓk =n

which already indicates a connection of the r-Lah numbers with the random 2r-compositions.

6.3. Definition of the (r, s)-Lah distribution.

Definition 6.1. An (r, s)-Lah distribution Lah[n, k]r, s with parameters n ∈ N, k ∈ {0, . . . , n} and
r, s ∈ [0, ∞) (where the case k = r = s = 0 is always excluded) is a discrete probability measure
defined by
   
1 n j
Lah[n, k]r, s ({j}) = , j ∈ {k, k + 1, . . . , n}. (6.6)
L(n, k) r+s j r k s
2

Definition 6.2. We call (n, k, r, s) an admissible quadruple of parameters for the (r, s)-Lah distri-
bution if n ∈ N, k ∈ {0, . . . , n}, r, s ∈ [0, ∞) and max{k, r, s} > 0.

We exclude the case k = r = s = 0 in which L(n, 0)0 = 0, n ∈ N, since (−1)! = ∞. Equation (6.6)
indeed defines a probability distribution due to identity (6.5).
Now we review several particular cases of the (r, s)-Lah distribution that have already appeared
in the literature.

Example 6.3. If k = 0 and s ≥ 0, then in view of 0j s = sj , the distribution takes the form


 
1 n j
Lah[n, 0]r, s ({j}) = n↑
s , j ∈ {0, 1, . . . , n}.
(r + s) j r
Thus, Lah[n, 0]r, s is the r-Stirling distribution of the first kind r-Stir1[n; s]. In particular, if also
s = 0, then Lah[n, 0]r, 0 is the degenerate at 0 distribution for each r > 0.
32
Example 6.4. Let now k = 1 and s = 0. Note that, for j ∈ {1, . . . , n},
       
(6.2) 1 n j (1.1) 1 n n+r 1 n
Lah[n, 1]r, 0 ({j}) = n+r−1

= n+r−1

= n↑
.
n! r j r 1 0 n! r j r n (r + 1) j r
Thus, if Z ∼ r-Stir1[n; 1], then Lah[n, 1]r, 0 is the conditional distribution of Z given {Z > 0}.
Example 6.5. Motivated by [39, 40] the (r, s)-Lah distribution with r = s = 0 was investigated
in [36]. This distribution has a combinatorial interpretation in terms of fragmented permutations.
A fragmented permutation of n elements with k fragments is a partition of the set {1, . . . , n}
into k nonempty blocks together with the structure of permutation on each block. The number of
fragmented permutations with k blocks is the Lah number L(n, k)0 . Consider a random fragmented
permutation with k blocks sampled uniformly from the set of all such fragmented permutations.
Then, the number of cycles of this fragmented permutation has the (r, s)-Lah distribution with
parameters r = s = 0.
Example 6.6. Consider a complete graph on n + 1 nodes in which the edges have i.i.d. unit
exponential weights. The union of shortest paths from some node of this graph to the k uniformly
chosen nodes is a random tree. The number of edges in this tree, denoted by Hn+1 (k), is analyzed
in [24]. According to Theorem 2.1 of that paper
   
k! n + 1 j
P[Hn+1 (k) = j] = , j ∈ {k, k + 1, . . . , n}.
n! nk j + 1 0 k 0

Note that nj = n+1

   
1 j+1 0 . This means that Hn+1 (k) has the same distribution as Lah[n, k]1, 0 . We
will return to this example in Section 6.5.
6.4. Generating polynomial of the (r, s)-Lah distribution and a recursion for the prob-
ability mass function. The next result identifies the generating polynomial of the (r, s)-Lah
distribution (up to a normalizing factor). It is a generalization of [38, Lemma 3.3 or Lemma 3.6].
Lemma 6.7. For any n ∈ N, k ∈ {0, . . . , n}, r, s ≥ 0 and t ∈ C,
n    
X n j n! 
k 
tj = [xn ] (1 − x)−t − 1 (1 − x)−(st+r) .
j r k s k!
j=k

Proof. We compute the Taylor series using the exponential generating functions (1.3) and (1.4):
  
1 −t

k −(st+r) −r 1 −t log(1−x) k
−st log(1−x)
(1 − x) − 1 (1 − x) = (1 − x) e −1 e
k! k!
∞
j ! ∞  
!
−t log(1 − x) (− log(1 − x))j
 
(1.4) −r
X j X j j −r
= (1 − x) = t (1 − x)
k s j! k s j!
j=k j=k
 
∞   ∞
X n xn
  ∞ ∞
X X 1 n  j 

(1.3) X j j
= t  = tj xn .
k s j r n! n! j r k s
j=k n=j j=k n=j

Taking the coefficient of xn gives the claimed identity. □

Corollary 6.8. Fix some n ∈ N and k ∈ {0, . . . , n}. Let r, s ∈ N0 and m ∈ N0 satisfy m <
(n + r)/(k + s). Then,
n    
X n j
(−m)j = 0.
j r k s
j=k
33
Proof. We apply Lemma 6.7 with t := −m. The degree of the polynomial (1 − x)m − 1)k (1 − x)sm−r
is km + sm − r < n under our assumption on m. Hence, the coefficient of xn vanishes. □
Corollary 6.9. Fix some n ∈ N, k ∈ {0, . . . , n} and r ∈ R. Let m ∈ N0 satisfy m < n/k. Then,
n    
X n j
(−m)j = 0.
j r k r
j=k

Proof. If r ∈ Nis sufficiently large, then m < (n + r)/(k + r) and we can apply Corollary 6.8 with
n
j
r = s. Since j and k r are polynomials in r, the identity extends to all real r. □
r
Now we derive recurrence relations for the probability mass function and the generating poly-
nomial of Lah[n, k]r, s . The triangles of generalized Stirling numbers satisfy the recursions [54,
p. 1661]:
           
n n−1 n−1 j j−1 j−1
= (n + r − 1) + , = + (k + s) .
j r j r j−1 r k s k−1 s k s
Upon multiplying these identities we conclude that
               
n j n−1 j n−1 j−1 n−1 j−1
= (n + r − 1) + + (k + s) . (6.7)
j r k s j r k s j−1 r k−1 s j−1 r k s
After simple manipulations with binomial coefficients we obtain from (6.7) the following trivariate
recursion for the probability mass function of the (r, s)-Lah distribution. Put
pn,k (j) = P[Lah[n, k]r, s = j].
Proposition 6.10. For j ∈ {k, k + 1, . . . , n}, the following trivariate recursive formula holds true
(n + r − 1)(n − k) (k + s)(n − k)
pn,k (j) = pn−1,k (j) + pn−1,k (j − 1)
n(n + r + s − 1) n(n + r + s − 1)
k(k + r + s − 1)
+ pn−1,k−1 (j − 1). (6.8)
n(n + r + s − 1)
Multiplying both sides of (6.8) by tj and summing over j yields a bivariate recursion for the
generating functions
n
X
Gn,k (t) := pn,k (j)tj = E[tLah[n, k]r, s ],
j=k
which is,
 
(n + r − 1)(n − k) (k + s)(n − k) k(k + r + s − 1)
Gn,k (t) = Gn−1,k (t) + t + tGn−1,k−1 (t).
n(n + r + s − 1) n(n + r + s − 1) n(n + r + s − 1)
(6.9)

6.5. Stochastic representation in terms of random recursive trees. It is known, see [16,
Theorem O1] or [62], that in RRT[n] both the degree of the root and the distance from the root to
node n are Stir1[n; 1]-distributed. The next theorem unifies and generalizes both claims.
Theorem 6.11. Consider a random tree Hoppe[n; r + s] with root 0 of weight r + s and nodes
1, . . . , n. Think of the root as being divided into two components with weights r and s which means
that if some node is attached to the root, it chooses one of these components with probabilities
r/(r + s) and s/(r + s). Let An,k be a random subset of {1, . . . , n} consisting of k nodes sampled
34
uniformly without replacement. Let Cn,k be the set of all nodes which are attached directly to the
weight s component of the root. If Tn,k denotes the minimal subtree rooted at 0 and containing all
nodes from An,k ∪ Cn,k , then the number of edges in Tn,k is distributed according to Lah[n, k]r, s .
Before proving the theorem we consider some special cases.
Corollary 6.12. The minimal subtree of a random tree Hoppe[n; r] spanned by a random uniform
subset of k nodes (and the root) has Lah[n, k]r, 0 edges.
Proof. Choose s = 0 in Theorem 6.11 (meaning that the set Cn,k is empty). □
In the case r = 1 and s = 0, Theorem 6.11 recovers Proposition 3.1 in [24]. Note that in [24] the
results are stated in terms of the so-called shortest path trees; these can be reduced to the RRT’s,
as explained in [67, Sections 16.2, 16.3], [66, Sections 1,2] and [68, Section 6.2]; see also [24, 69].
In the case where r > 0 is arbitrary and k = 1, recalling a characterization of Lah[n, 1]r, 0 given
in Example 6.4 we arrive at the following
Corollary 6.13. The distance from the root of a random tree Hoppe[n; r] to a node Vn picked
uniformly at random from the set of nodes {1, . . . , n} has the same distribution as Z ∼ r-Stir1[n; 1]
conditioned on Z ̸= 0.
In the next corollary we provide a simple representation of the standard (r = s = 0) Lah
distribution.
Corollary 6.14. Consider a random recursive tree with root denoted for the time being by 1 and
nodes 2, . . . , n. Then, the number of nodes in a subtree generated by the root and a uniform
random subset of k nodes from {1, . . . , n} (note that 1 may be in the subset) is distributed according
to Lah[n, k]0, 0 .
Proof. We consider a Hoppe[n; r]-tree with r tending to zero and apply Corollary 6.12. Note that
node 1 is always attached to 0, while nodes 2, . . . , n are attached to 0 with probability tending
to zero as r → 0. Hence, if we remove 0 and the edge between 0 and 1, then with probability
tending to one we obtain the random recursive tree described in Corollary 6.14. It remains to
apply Corollary 6.12 and observe that the number of nodes in the RRT is the number of edges plus
the removed edge between 0 and 1. □
In the case k = 1, Corollary 6.14 recovers a known result on the insertion depth in the RRT;
see [62, Theorem 2] or [65].
Example 6.15. The case k = 0 of Theorem 6.11 is straightforward. In this scenario An,k is empty
and Cn,k has the mixed binomial distribution Bin[Stir1[n; r + s]; s/(r + s)], where Stir1[n; r + s] is
the degree distribution of the root in Hoppe[n; r+s]. By part (ii) of Proposition 3.7 Bin[Stir1[n; r+
s]; s/(r + s)] = r-Stir1[n; s] and by Example 6.3 r-Stir1[n; s] = Lah[n, 0]r, s in full accordance with
the k = 0 case of Theorem 6.11.
Proof of Theorem 6.11. Replacing n by n + 1, we can write recurrence relation (6.9) for Gn,k (t) =
E[tLah[n, k]r, s ] in the form
  
k n+r st
Gn+1,k (t) = 1 − Gn,k (t) + Gn,k (t) (6.10)
n+1 n+r+s n+r+s
 
k k−1+r+s n−k+1
+t· Gn,k−1 (t) + Gn,k (t) . (6.11)
n+1 n+r+s n+r+s
35
We show that the generating function of the number of edges in Tn,k satisfies the same recurrence
relation as Gn,k . Our argument generalizes the one in [24]. Consider the tree Hoppe[n + 1; r + s].
With probability k/(n + 1), node n + 1 belongs to An+1,k . Denote this event by Q. Given Q, there
two possibilities:
Case 1. If the node n + 1 is attached to a node from the set An,k−1 or to the root with weight r + s,
then the tree Tn+1,k differs from Tn,k−1 by one edge containing the node n + 1. This additional
edge accounts for a factor t in (6.11). The conditional (given Q) probability of the event occurring
in Case 1 is (k − 1 + r + s)/(n + r + s).
Case 2. If the node n + 1 is attached to some node v from the set {1, . . . , n} \ An,k , then the tree
Tn+1,k can be constructed as follows. Construct the minimal tree containing the root and the set
An,k−1 ∪ {v} ∪ Cn,k and add to this tree the edge connecting n + 1 and v. Then, Tn+1,k differs by
one edge from the tree having the same law as Tn,k , which again gives a factor t. The conditional
(given Q) probability of the event occurring in Case 2 is (n − k + 1)/(n + r + s). Altogether, these
considerations explain the term (6.11).
Now, it is possible that the node n+1 does not belong to An+1,k , which happens with probability
P[Qc ] = 1 − k/(n + 1). Given Qc , there are again two possibilities:
Case 3. If the node n + 1 is attached to the weight s component of the root, then it belongs to
Cn,k and the tree Tn+1,k differs from Tn,k by one edge. The conditional probability (given Qc ) of
this event is s/(n + r + s).
Case 4: If the node n + 1 is attached to any node in {1, . . . , n} or to the weight r component of the
root, then Tn+1,k = Tn,k . The conditional probability (given Qc ) of this event is (n + r)/(n + r + s).
Altogether, Cases 3 and 4 explain (6.10).
It remains to note that on the boundary k = 0 the generating function of the number of edges
in Tn,0 also coincides with Gn,0 (t) = E[tLah[n, 0]r, s ] as we have seen in Example 6.15. □

6.6. Stochastic representation in terms of random compositions. The next proposition pro-
vides a representation of the (r, s)-Lah distribution as a sum of conditionally independent random
variables. It will be the starting point in our asymptotic analysis of the (r, s)-Lah distribution.
(n,k) (n,k) (n,k)
Proposition 6.16. Let (b0 , b1 , . . . , bk ) be a random (r + s)-composition of n. Further,
(n,0) (n,1) (n,k) (n,k)
let (Zr,s , Zr,s , . . . , Zr,s ) be a random vector such that, conditionally on the event {b0 =
(n,k) (n,k)
b0 , b1 = b1 , . . . , bk = bk }, where (b0 , b1 , . . . , bk ) is an arbitrary fixed incomplete composition of
n, the random variables
(n,0) (n,1) (n,k)
Zr,s ∼ r-Stir1[b0 ; s], Zr,s ∼ Stir1[b1 ; 1], ..., Zr,s ∼ Stir1[bk ; 1]
are independent, that is,
 
h i bj
(n,k) (n,k) (n,k)
X
(n,j)
P Zr,s = ℓ b0 = b0 , b1 = b1 , . . . , bk = bk = P  Bern[1/m] = ℓ
m=1

for all j = 1, . . . , k and

" b #
h i 0
(n,k) (n,k) (n,k)
X
(n,0)
P Zr,s = ℓ b0 = b0 , b1 = b1 , . . . , bk = bk = P Bern[s/(r + s + m − 1)] = ℓ .
m=1
36
Here, all Bernoulli random variables are assumed independent. Then,
(n,0) (n,1) (n,k)
Zr,s + Zr,s + · · · + Zr,s ∼ Lah[n, k]r, s .
Example 6.17. Taking k = 0 gives Lah[n, 0]r, s = r-Stir1[n; s].
Example 6.18. For r = s = 0, Proposition 6.16 reduces to the representation of the Lah distribu-
tion given in [36, Proposition 2.3].
Example 6.19 (Cumulative degree of roots in the multinomial Hoppe tree). Consider a random
forest MultiHoppe[n − k; θ; ((r + s)/θ, 1/θ, . . . , 1/θ)] with root0 having weight r + s and k roots
root1 , . . . , rootk having weight 1 each. The total weight of all roots is θ = k + r + s. According
to Proposition 5.7 and Proposition 2.8, the vector encoding the sizes of connected components
(not counting root0 but counting root1 , . . . , rootk ) has the same distribution as a random (r + s)-
composition of n. Let Dn;j be the degree of rootj , for j = 0, . . . , k. Then, Proposition 6.16 implies
that
Dn;1 + · · · + Dn;k + Bin[Dn;0 ; s/(r + s)] ∼ Lah[n, k]r, s .
For r = 0 the right-hand side simplifies to Dn;0 + Dn;1 + · · · + Dn,k .
(n,j)
Example 6.20. For r = s = 1/2, the conditional distribution (as in Proposition 6.16) of Zr,s
Pbj
with j ∈ {1, . . . , k} is the same as m=1 Bern[1/m], which is the distribution of the number of
records in bj i.i.d. observations from a continuous distribution, while the conditional distribution of
(n,0)
Zr,s is the same as bm=1
P0
Bern[1/(2m)]. The distribution of the random 1-composition is uniform
on the set of all incomplete compositions of n. Motivated by the analysis of positive hulls of random
walks carried out in [21] (see also [20]) the case r = s has already been investigated in [38].
Proof of Proposition 6.16. Recall from Lemma 6.7 and (6.2) that
1 
k 
Gn,k (t) = E[tLah[n, k]r, s ] = n+r+s−1
[xn ] (1 − x)−t − 1 (1 − x)−(st+r) . (6.12)
k+r+s−1
(n,0) (n,1) (n,k)
From the definition of the random variables Zr,s , Zr,s , . . . , Zr,s it follows that
h (n,j) i tbj ↑
(n,k) (n,k) (n,k)
E tZr,s b0 = b0 , b1 = b1 , . . . , bk = bk = , j ∈ {1, . . . , k}
bj !
and
h (n,0) i (st + r)b0 ↑
(n,k) (n,k) (n,k)
E tZr,s b0 = b0 , b1 = b1 , . . . , bk = bk = .
(r + s)b0 ↑
By conditional independence, this gives
h (n,0) (n,1) (n,k)
i tb1 ↑ tbk ↑ (st + r)b0 ↑
(n,k) (n,k) (n,k)
E tZr,s +Zr,s +···+Zr,s b0 = b0 , b1 = b1 , . . . , bk = bk = ·· · ·· · . (6.13)
b1 ! bk ! (r + s)b0 ↑
Also, by the definition of the random (r + s)-composition (Definition 5.2),
(n,k) (n,k) (n,k) (r + s)b0 ↑ /b0 !
P[b0 = b0 , b1 = b1 , . . . , bk = bk ] = n+r+s−1

. (6.14)
k+r+s−1
Using the formulas
∞ ∞
−t
X xb −(st+r)
X xb0
(1 − x) −1= b↑
t and (1 − x) = (st + r)b0 ↑ ,
b! b0 !
b=1 b0 =0
37
we obtain
1 X tb1 ↑ tbk ↑ (st + r)b0 ↑
Gn,k (t) = n+r+s−1

...
k+r+s−1
b1 ! bk ! b0 !
(b0 ,b1 ,...,bk )
X tb1 ↑ tbk ↑ (st + r)b0 ↑ (r + s)b0 ↑ /b0 !
= · ··· · · n+r+s−1
,
b1 ! bk ! (r + s)b0 ↑ k+r+s−1
(b0 ,b1 ,...,bk )

where both sums are taken over all incomplete compositions of n. Recognizing on the right-hand side
(n,0) (n,1) (n,k) (n,k) (n,k) (n,k)
the conditional expectation of tZr,s +Zr,s +···+Zr,s given {b0 = b0 , b1 = b1 , . . . , bk = bk },
see (6.13), and the probability of this event, see (6.14), we apply the total expectation formula to
obtain
(n,0) (n,1) (n,k)
Gn,k (t) = E[tLah[n, k]r, s ] = E[tZr,s +Zr,s +···+Zr,s ].
The proof is complete. □

Using Corollary 5.12 and recalling the representation of the standard (r = s = 0) Lah distribution
given in [36, Proposition 2.3] we obtain the following.
(n,k) (n,k) (n,k)
Corollary 6.21. Let (b0 , b1 , . . . , bk ) be a random (r+s)-composition and (Lah[m, k]0, 0 )m≥k
a sequence of random variables having the standard Lah distributions, which is independent of
(n,0) (n,k) (n,0)
(Zr,s , b0 ), where Zr,s is as defined in Proposition 6.16. Then,
(n,0) (n,k)
Lah[n, k]r, s has the same distribution as Zr,s + Lah[n − b0 , k]0, 0 .

Corollary 6.21 is designed to derive limit theorems for Lah[n, k]r, s from the known limit theorems
for Lah[n, k]0, 0 .

6.7. Expectation and variance. We employ the notation

L(n, k, r, s) := E[Lah[n, k]r, s ].
The following explicit formula can be found in Theorem 3.1 of [38]
nk
L(n, k, 0, 0) = E[Lah[n, k]0, 0 ] = (Hn − Hk−1 ), (6.15)
n − (k − 1)
(r,s)
where Hn := 1 + 1/2 + · · · + 1/n is the n-th harmonic number. For r, s ≥ 0, put H0 := 0 and,
(r,s) Pn −1 (r,0)
for n ∈ N, Hn := m=1 s(r + s + m − 1) . Plainly, Hn = 0 for all n ∈ N0 . Corollary 6.21
immediately yields
(r,s) (n,k)
L(n, k, r, s) = E[H (n,k) ] + E[L(n − b0 , k, 0, 0)],
b0
(n,k) (n,k)
where (b0 , . . . , bk ) is a random (r + s)-composition. However, this formula is of little use for
practical purposes. Note that Proposition 6.16 implies that
(r,s)
L(n, k, r, s) = E[H (n,k) ] + kE[Hb(n,k) ].
b0 j

We do not derive exact formulas for either L(n, k, r, s) or the variance, although it can be done
relatively straightforwardly by generalizing the proofs given in [38, Theorems 3.1, 3.7]. Instead,
our focus is on asymptotic estimates, which are crucial for proving limit theorems for the Lah
distributions in subsequent sections.
38
Lemma 6.22. For r, s ≥ 0, the following estimate holds true
 r + s + (n − 1)  s
+
Hn(r,s) − s log ≤ , n ∈ N0 . (6.16)
r+s r+s
Proof. For n ∈ N0 ,
n Z (n−1)+  r + s + (n − 1) 
X s s s dx s +
Hn(r,s) = ≤ + = + s log .
r+s+m−1 r+s 0 r + s + x r + s r + s
m=1
This together with an analogous estimate from below completes the proof of (6.16). □
(n,k) (n,k)
Lemma 6.23. Fix some r, s ≥ 0 and let (b0 , . . . , bk ) be a random (r+s)-composition. Assume
that n, k → ∞ such that k/n → 0. Then, for each fixed a ≥ 1,
h k k
(r,s)
a
X
a i hX
a i
E H (n,k) + Hb(n,k) ≃ E Hb(n,k) ≃ k(log(n/k))a , n, k → ∞. (6.17)
b0 j j
j=1 j=1

In particular,
k
hX i
L(n, k, r, s) ≃ E Hb(n,k) ≃ k log(n/k), n → ∞. (6.18)
j
j=1
Furthermore,
k
hX i
Var Hb(n,k) = O(k), n→∞ (6.19)
j
j=1
and
h k i
(r,s)
X
Var H (n,k) + Hb(n,k) = O(k), n → ∞. (6.20)
b0 j
j=1
(r,s) (r,0)
Proof. When dealing with H (n,k) we always assume that s > 0, for H (n,k) = 0 a.s.
b0 b0
(r,s)
For a proof of (6.17) it suffices to check that E[(Hb(n,k) )a ] ≃ (log(n/k))a and E[(H (n,k) )a ] =
1 b0
O(log(n/k)) as n, k → ∞. In view of (5.11),
(n,k)
log b1 P
→ 1, n, k → ∞.
log(n/k)
This in combination with |Hn − log n| ≤ C for all n ∈ N and a constant C > 0 entails
(Hb(n,k) )a
1 P
→ 1, n, k → ∞
(log(n/k))a
and thereupon by Fatou’s lemma lim inf n,k→∞ (E[(Hb(n,k) )a ]/(log(n/k))a ) ≥ 1. To prove the con-
1
verse inequality for the upper limit, note that x 7→ (log x + C)a is concave on [1, ∞) for each
C > a − 1. In view of this we take C sufficiently large and write with the help of Jensen’s inequality
(n,k) (n,k)
E[(Hb(n,k) )a ] ≤ E[(log b1 + C)a ] ≤ (log E[b1 ] + C)a
1
(n,k)
= (log E[(n − b0 )/k] + C)a ≤ (log(n/k) + C)a .
Using (6.16), we conclude that
 a  r + s + (n − 1) 
s +
(Hn(r,s) )a ≤2 a−1
+ s log a a
.
r+s r+s
39
Further, Jensen’s inequality and formula (5.10) with r replaced by r + s yields
" #

s
a   r + s + (b(n,k) − 1)  a
(r,s) +
E[(H (n,k) )a ] − 2a−1 ≤ 2a−1 sa E log 0
+C
b0 r+s r+s
" #
  r + s + E[b(n,k) ]  a
a−1 a 0
≤2 s log +C
r+s
 a 
a−1 a
n−k 
=2 s log 1 + +C . (6.21)
k+r+s

The right-hand side is O((log(n/k))a ). The proof of (6.17) is complete.

(r,s)
We now turn to the proof of (6.19) and (6.20). The functions x 7→ H⌊x⌋ and x 7→ H⌊x⌋ are
nondecreasing on [0, ∞) and [1, ∞), respectively. An application of Corollary 5.10, with r replaced
by r + s, then yields, for j ∈ {1, . . . , k},
 (r,s)   (r,s)   
E H (n,k) Hb(n,k) ≤ E H (n,k) E Hb(n,k)
b0 j b0 j

and, for i, j ∈ {1, . . . , k}, i ̸= j,

     
E Hb(n,k) Hb(n,k) ≤ E Hb(n,k) E Hb(n,k) .
i j i j

Thus, both (6.19) and (6.20) follow provided we can prove

(r,s)
Var H (n,k) = O(1) and Var Hb(n,k) = O(1), n, k → ∞. (6.22)
b0 1

To check this, we show that, for each p ≥ 1,

(n,k)
p
sup E log r + s + (b0 − 1)+ − log(n/k) < ∞ and
n,k∈N,k≤n
(n,k) p
sup E log b1 − log(n/k) < ∞. (6.23)
n,k∈N,k≤n

Our argument is a slight adaptation of the proof of formula (5.33) in [36]. We start by proving the
second inequality in (6.23). In view of |x| = x+ + x− for x ∈ R and (x + y)p ≤ 2p−1 (xp + y p ) for
x, y ≥ 0 it suffices to show that

(n,k) p (n,k) p
sup E[log+ ((k/n)b1 )] < ∞ and sup E[log− ((k/n)b1 )] < ∞.
n,k∈N,k≤n n,k∈N,k≤n

(n,k)
The former inequality follows from [log+ x]p = O(x) as x → ∞ and the estimate (k/n)E[b1 ] =
(n,k)
(n − E[b0 ])/n ≤ 1. To prove the latter inequality, we first conclude with the help of (5.7) that
 (n,k) 
P b1 =j k+r+s−2
 (n,k)  =1+ >1
P b1 =j+1 n−j−k+1
40
  (n,k) 
for all k > (2 − r − s)+ and all j ∈ {1, . . . , n − k}. Thus, the sequence j 7→ P b1 = j is strictly
decreasing. This entails
⌊n/k⌋ ⌊n/k⌋
(n,k) (n,k) (n,k)
X X
E[log− ((k/n)b1 )]p = [log− (jk/n)] p
P[b1 = j] ≤ P[b1 = 1] [log− (jk/n)]p
j=1 j=1
⌊n/k⌋ Z 1
k+q−1 X
= [log− (jk/n)]p → | log x|p dx < ∞, n, k → ∞.
n+q−1 0
j=1

The proof of the second inequality in (6.23) is complete.

Turning to the first inequality in (6.23), we conclude that it is enough to show that
(n,k)
sup E[log+ ((k/n)(r + s + (b0 − 1)+ )]p < ∞ and
n,k∈N,k≤n
(n,k)
sup E[log− ((k/n)(r + s + (b0 − 1)+ ))]p < ∞.
n,k∈N,k≤n

The first of these is secured by

(n,k) (n,k)
(k/n)E[r + s + (b0 − 1)+ ] ≤ r + s + (k/n)E[b0 ]
= r + s + (k/n)((r + s)(n − k)/(k + r + s)) ≤ 4(r + s).
We have used (5.10) for the equality. As for the second, we write
(n,k) (n,k)
E[log− ((k/n)(r + s + (b0 − 1)+ ))]p = [log− ((r + s)(k/n))]p P[(b0 − 1)+ = 0]
⌊n/k−(r+s)⌋
(n,k)
X
+ [log− ((k/n)(r + s + j))]p P[b0 = j + 1] =: I1 (n, k) + I2 (n, k).
j=1

(n,k)
In view of (5.6), P[(b0 − 1)+ = 0] ≃ (r + s + 1)(k/n)r+s as n, k → ∞. This implies that the first
summand I1 vanishes. Now we analyze I2 . Since
(n − j − 1)! (n − j − 1 − k)! n−j−1
= >1
(n − j − k)! (n − j − 2)! n−j−k
n−j−1


for all k ≥ 2 and all j, n ∈ N satisfying n − j − k ≥ 1, we conclude that the sequence j 7→ k−1
is decreasing. Thus, for j ≥ 1,
n−j−2 n−3



 k r+s
k−1 k−1
n+r+s−1

≤ n+r+s−1

≃ , n, k → ∞.
k+r+s−1 k+r+s−1
n
By a standard estimate for the gamma function, there exists a constant cr,s > 0 such that, for all
j ≥ 1,
 
r+s+j
≤ cr,s j r+s−1 .
1+j
Combining fragments together and recalling (5.6) we infer
 k  ⌊n/k−(r+s)⌋
X  kj r+s−1
I2 (n, k) = O [log− ((k/n)(r + s + j))]p = O(1), n, k → ∞.
n n
j=1
41
The last equality is justified by a Riemann approximation
⌊n/k−(r+s)⌋ Z 1
k X
p
 kj r+s−1
lim [log− ((k/n)(r + s + j))] = | log x|p xr+s−1 dx < ∞.
n,k→∞ n n 0
j=1

This completes the proof of (6.23).

We only show how to obtain the first relation in (6.22), the argument for the second is analogous.
In view of (6.16) and the first inequality in (6.23),
(r,s)  (r,s) 2
E[H (n,k) ] − s log(n/k) = O(1) and E H (n,k) − s log(n/k) = O(1), n, k → ∞.
b0 b0

The latter entails

 (r,s) 2 (r,s)
E H (n,k) = 2s log(n/k)E[H (n,k) ] − (s log(n/k))2 + O(1), n, k → ∞,
b0 b0

whence
(r,s)  (r,s) 2  (r,s) 
2 (r,s)
Var [H (n,k) ] = E H (n,k) − E H (n,k) = −(EH (n,k) − s log(n/k))2 + O(1) = O(1), n, k → ∞.
b0 b0 b0 b0

The proof is complete. □

Remark 6.24. A perusal of the proof of (6.17) reveals that, for fixed k,
h i
E Hb(n,k) = s log n + O(1), n → ∞,
0

whence
L(n, k, r, s) = (k + s) log n + O(1), n → ∞.

7. Limit theorems for the (r, s)-Lah distribution

In this concluding section we prove central limit theorems for the (r, s)-Lah distribution. Not only
do these theorems generalize the special cases r = s = 0 and r = s investigated  in [36] and [38], but
2

they also treat asymptotic regimes not covered in those papers. Let N 0; σ denote a normal ran-
(n,k) (n,k) (n,k)
dom variable with zero mean and variance σ 2 > 0. Throughout this section (b0 , b1 , . . . , bk )
denotes the random (r + s)-composition.

Theorem 7.1 (CLT in the constant k regime). Let k ∈ N0 and r, s ∈ [0, ∞) be fixed and such that
max{k, s} > 0. Then,
Lah[n, k]r, s − (k + s) log n d
p −→ N [0; 1] .
(k + s) log n n→∞

Proof. We assume that r + s > 0. The proof hinges on Corollary 6.21 and a known CLT in the
case of constant k and r = s = 0, see Theorem 4.2 in [36]. A classic functional limit theorem for
independent Bernoulli variables together with Theorem 4.2 in [36] imply
 
P⌊nt ⌋ !
(n,k) (n,k)
m=1 Bern[s/(r√+ s + m − 1)] − ts log n Lah[n − b0 , k]0, 0 − k log(n − b0 ) 
 , q
s log n (n,k)
t≥0 k log(n − b0 )
d
−→ ((B(t))t≥0 , N [0; 1]) (7.1)
n→∞
42
with the components on both sides being independent and (B(t))t≥0 being a standard Brownian
(n,k) d
motion. According to Proposition 5.15, b0 /n −→ Beta [r + s; k], whence
n→∞

(n,k) (n,k)
log+ b0 P log(n − b0 ) − log n P
→ 1 and √ → 0, n → ∞. (7.2)
log n log n
Using continuity of the composition operation and Slutsky’s lemma we conclude from (7.1) and (7.2)
that
(n,0) (n,k)
!
(n,k) √ √
Zr,s − s log+ b0 Lah[n − b0 , k]0, 0 − k log n d

√ , √ −→ sB(1), kN [0; 1] . (7.3)
log n log n n→∞

By adding the coordinates and using that, see Proposition 5.15,

(n,k)
log+ b0 − log n P
√ → 0, n → ∞,
log n
we obtain the desired statement. The proof is complete. □

Theorem 7.1 was previously known in case r = s > 0 from [38] and in case r = s = 0 from [36].
The next theorem is new even in the standard scenario r = s = 0.

Theorem 7.2 (CLT in the intermediate regime). Fix some r, s ≥ 0 and let n → ∞, k → ∞ such
that k/n → 0. Then,
Lah[n, k]r, s − L(n, k, r, s) d
p −→ N [0; 1] .
k log(n/k) n,k→∞

Proof of Theorem 7.2. According to Proposition 6.16, it is sufficient to show that

(n,0) (n,k)
Zr,s + · · · + Zr,s − L(n, k, r, s) d
p −→ N [0; 1] .
k log(n/k) n,k→∞

(n,k) (n,k) (n,k)

For each n ∈ N and k ∈ N, denote by Gn,k the σ-algebra generated by b0 , b1 , . . . , bk . We
first prove that, conditionally on Gn,k ,
Pk (n,j)  (n,j) 

j=0 Zr,s − E Zr,s Gn,k d
q −→ N [0; 1] . (7.4)
Pk (n,j) n,k→∞
Var [ j=0 Zr,s |Gn,k ]

Conditionally on Gn,k , the left-hand side is a sum of independent centered random variables with
finite third absolute moments. In view of this, we will prove the last limit relation by an application
of the Berry–Esseen inequality:
P 
k (n,j)  (n,j) 

j=0 Z r,s − E Z r,s G n,k
In,k := sup P  q ≤ x − P[N [0; 1] ≤ x]
x∈R
P k (n,j)
Var [ j=0 Zr,s |Gn,k ]
Pk (n,j)  (n,j)
Gn,k |3 |Gn,k ]

j=0 E[|Zr,s − E Zr,s
≤A (n,j)
(Var [ kj=0 Zr,s |Gn,k ])3/2
P

43
for a deterministic constant A > 0 which does not depend on n, nor k. To complete the proof
of (7.4) it is enough to show that
Pk (n,j)  (n,j)
Gn,k |3 |Gn,k ] P

j=0 E[|Zr,s − E Zr,s
(n,j)
→ 0, n, k → ∞. (7.5)
(Var [ kj=0 Zr,s |Gn,k ])3/2
P

As a preparation, we prove that

 Pk (n,j) 
Var j=0 Zr,s Gn,k P
→ 1, n, k → ∞. (7.6)
k log(n/k)
(n,0)
In what follows we assume that s > 0, for otherwise Zr,s = 0 a.s. In view of (6.16),
(n,k)
b0
 (n,0)  X s  s 
(r,s)
Var Zr,s Gn,k = 1− ≤ H (n,k)
r+s+m−1 r+s+m−1 b0
m=1
 r + s + (b(n,k) − 1)  s
0 +
≤ s log + .
r+s r+s
(n,k) P
As a consequence of (5.11), log(b0 − 1)+ / log(n/k) → 1 as n, k → ∞. This shows that
 (n,0) 
Var Zr,s Gn,k P
→ 0, n, k → ∞.
k log(n/k)
Since, for j ∈ N,
(n,k)
bj
X 1 1
(n,j)
∈ [Hb(n,k) − π 2 /6, Hb(n,k) ],
 
Var Zr,s Gn,k = 1−
m m j j
m=1
relation (7.6) is equivalent to
Pk
j=1 Hb(n,k) P
j
→ 1, n, k → ∞.
k log(n/k)
The latter is secured by Lemma 6.23 and Chebyshev’s inequality. Thus, for the proof of (7.5) it
remains to check
Pk (n,j)  (n,j)
Gn,k |3 |Gn,k ] P

j=0 E[|Zr,s − E Zr,s
→ 0, n, k → ∞. (7.7)
(k log(n/k))3/2
Invoking once again Proposition 6.16 and using Rosenthal’s inequality, see Theorem 9.1 on p. 152
in [23], we obtain for an absolute constant B > 0 which does not depend on n and k
(n,k)
 bX
j
2 3/2
3
E Bern[m−1 ] − m−1
 (n,j)  (n,j)   
E Zr,s − E Zr,s Gn,k Gn,k ≤ B
m=1
(n,k)
bj
X 3 
3/2
E Bern[m−1 ] − m−1 |


+ ≤ B Hb(n,k) + Hb(n,k) a.s.
j j
m=1
for j ≥ 1 and analogously
 (n,0)  (n,0)  3  (r,s)
3/2 (r,s)

E Zr,s − E Zr,s Gn,k Gn,k ≤ B H (n,k) + H (n,k) a.s.
b0 b0
44
Therefore, (7.7) is secured by Eq. (6.17) with a = 3/2 and Markov’s inequality.
We have proved that (7.4) holds conditionally on Gn,k . By Lebesgue’s dominated convergence
theorem, it also holds unconditionally. Furthermore, in view of (7.6),
Pk (n,j)  (n,j) 

j=0 Zr,s − E Zr,s Gn,k d
p −→ N [0; 1] .
k log(n/k) n,k→∞

Observe that
h k i k k
(r,s) (r,s)
X X  (n,j)  X
L(n, k, r, s) = E H (n,k) + Hb(n,k) , E Zr,s Gn,k = H (n,k) + Hb(n,k) .
b0 j b0 j
j=1 j=0 j=1

The proof finishes by an appeal to Eq. (6.20) in Lemma 6.23, which ensures that
Pk  (n,j) 
j=0 E Zr,s Gn,k − L(n, k, r, s) P
p → 0, n, k → ∞.
k log(n/k)
The proof is complete. □
Theorem 7.3 (CLT in the central regime). Fix some r, s ≥ 0 and let n → ∞, k → ∞ such that
k/n → α for some α ∈ (0, 1). Then,
α(α + 1) log α α2 log2 α
  
Lah[n, k]r, s − L(n, k, 0, 0) d α
√ −→ N 0; − + + .
n n,k→∞ 1−α (1 − α)2 (1 − α)3
Proof. Our proof relies on Corollary 6.21 in conjunction with the fact that in the central regime
(n,k)
b0 converges in distribution to the negative binomial distribution; see Proposition 5.16.
We first check that
(n,k)
Lah[n, k]0, 0 − Lah[n − b0 , k]0, 0 P
√ → 0, n, k → ∞. (7.8)
n
In view of Proposition 2.4 in [36] (or by using the couplings from Section 5.1) and conditional
Markov’s inequality it is sufficient to show that
(n,k)
L(n, k, 0, 0) − L(n − b0 , k, 0, 0) P
√ → 0, n, k → ∞. (7.9)
n
(n,k)
Indeed, (7.9) implies that (7.8) holds conditionally on b0 and, by the Lebesgue dominated con-
vergence theorem, also unconditionally.
Recall that L(n, k, 0, 0) is given explicitly in Eq. (6.15). Convergence (7.9) follows immediately
from the estimates, with a ∈ N0 fixed,
k
0 ≤ L(n, k, 0, 0) − L(n − a, k, 0, 0) ≤ (n(Hn − Hk−1 ) − (n − a)(Hn−a − Hk−1 ))
n−k+1
k
≤ (2aHn + n(Hn − Hn−a ))
n−k+1
  
k a
≤ 2aHn + n min Hn , ,
n−k+1 n−a
(n,k)
and the fact that b0 converges in distribution. The latter fact also ensures that
(n,0)
Zr,s P
√ → 0, n → ∞,
n
45
see Proposition 6.16. Combining this with (7.8) and appealing to Corollary 6.21 shows that the
claim of the theorem is equivalent to

α(α + 1) log α α2 log2 α

  
Lah[n, k]0, 0 − L(n, k, 0, 0) d α
√ −→ N 0; − + + .
n n,k→∞ 1−α (1 − α)2 (1 − α)3
But this is secured by Theorem 5.1 in [36]. The proof is complete. □

The next result extends Theorem 7.3 to the border case α = 1 in which the limiting variance in
Theorem 7.3 vanishes.

Theorem 7.4. Fix some r, s ≥ 0 and let n → ∞, k → ∞ such that k/n → 1 and n − k → ∞.
Then,
Lah[n, k]r, s − L(n, k, r, s) d
√ −→ N [0; 1/4] .
n−k n,k→∞

Before delving into a proof we briefly explain the idea. In the regime k/n → 1, the number of
(n,k) (n,k) (n,k)
blocks in the incomplete composition (b0 , b1 , . . . , bk ) is of the same magnitude as n. As
we will see, this means that the variability of the number of blocks of sizes > 2 or of size 1 is
(n,0) (n,1) (n,k)
negligible in the sense that they only contribute to Lah[n, k]r, s ∼ Zr,s + Zr,s + · · · + Zr,s
through the expectations. The principal contribution to the variance of Lah[n, k]r, s comes from
the blocks of size 2, and the number of these blocks turns out to be sufficiently close to n − k in
(n,j)
probability. Independently for each block of size 2 the corresponding sum Zr,s takes values 1 or 2
with probability 1/2 resulting in the limiting variance 1/4. We formalize and summarize the above
heuristics in a lemma.
(n,k) (n,k) (n,k)
Lemma 7.5. For m ∈ N, m ≤ n − k + 1 and a random (r + s)-composition (b0 , b1 , . . . , bk ),
let ρm (n, k) be the number of its (usual) blocks of size m, that is,
k
(n,k)
X
ρm (n, k) = 1[bj = m].
j=1

Let n → ∞, k → ∞ such that k/n → 1 and n − k → ∞. Then

(a) Var [ρ1 (n, k)] ≃ n−1 (n − k)2 = o(n − k) as n, k → ∞.
P
(b) Var [ρ2 (n, k)] ≃ 4n−1 (n − k)2 = o(n − k) and ρ2 (n, k)/(n − k) → 1 as n, k → ∞.
(c) n−k+1 Hm E[ρm (n, k)] = O(n−1 (n − k)2 ) = o(n − k) as n, k → ∞.
P
hm=3
Pn−k+1 i2
(d) E m=3 Hm (ρm (n, k) − E[ρm (n, k)]) = O(n−1 (n − k)2 ) = o(n − k) as n, k → ∞.

Proof. According to (5.7), for m ∈ N,

n−m+r+s−1


(n,k) k+r+s−2
E[ρm (n, k)] = kP[b1 = m] = k n+r+s−1

.
k+r+s−1

According to (5.8), for 1 ≤ i < j ≤ k,

n+r+s−2m−1


(n,k) (n,k) k+r+s−3
P[bi = m, bj = m] = n+r+s−1

.
k+r+s−1
46
This yields
(n,k) (n,k)
X
Var [ρm (n, k)] = E[ρm (n, k)] + 2 P[bi = m, bj = m] − (E[ρm (n, k)])2
1≤i<j≤k
(n,k) (n,k)
= E[ρm (n, k)] + k(k − 1)P[b1 = m, b2 = m] − (E[ρm (n, k)])2
n−m+r+s−1

n+r+s−2m−1

n−m+r+s−1

!2
k+r+s−2 k+r+s−3 k+r+s−2
=k n+r+s−1

+ k(k − 1) n+r+s−1

− k n+r+s−1

.
k+r+s−1 k+r+s−1 k+r+s−1

We will make a repeated use of these formulas without further reference and in conjunction with
the identities
       
α α α−1 α α−β+1 α
= and = .
β β β−1 β β β−1
(a) Plugging m = 1 yields

k(k + r + s − 1) k(k − 1)(k + r + s − 1)(k + r + s − 2) k 2 (k + r + s − 1)2

Var [ρ1 (n, k)] = + −
n+r+s−1 (n + r + s − 1)(n + r + s − 2) (n + r + s − 1)2
k(n − k)(n − k + r + s − 1)(k + r + s − 1) (n − k)2
= ≃ = o(n − k), n, k → ∞.
(n + r + s − 2)(n + r + s − 1)2 n

(b) Plugging m = 2 yields

(k + r + s − 1)(n − k)
Var [ρ2 (n, k)] = k
(n + r + s − 2)(n + r + s − 1)
(k + r + s − 1)(k + r + s − 2)(n − k − 1)(n − k)
+ k(k − 1)
(n + r + s − 4)(n + r + s − 3)(n + r + s − 2)(n + r + s − 1)
 2
(k + r + s − 1)(n − k)
− k
(n + r + s − 2)(n + r + s − 1)
≃ 4n−1 (n − k)2 = o(n − k), n, k → ∞.

Noting that
(k + r + s − 1)(n − k)
E[ρ2 (n, k)] = k ∼ n − k, n, k → ∞,
(n + r + s − 2)(n + r + s − 1)
the claim about the convergence in probability follows with the help of Chebyshev’s inequality.
(c) For m ≥ 3,
m−2
k(k + r + s − 1) Y n−k−i  n−k m−1
E[ρm (n, k)] = ≤k .
n+r+s−1 n+r+s−2−i n+r+s−2
i=0

It remains to note that

n−k+1
X  n−k m−1 (n − k)2
k Hm ∼ H3 , n, k → ∞,
n+r+s−2 n
m=3

by the Lebesgue dominated convergence theorem.

47
(d) By the triangle inequality in L2 , subadditivity of x 7→ x1/2 on [0, ∞) and the inequality (x+y)2 ≤
2x2 + 2y 2 for x, y ≥ 0
h n−k+1
X i2  n−k+1
X 2
E Hs (ρm (n, k) − E[ρm (n, k)]) ≤ Hm (Var [ρm (n, k)])1/2
m=3 m=3
 n−k+1 2
(n,k) (n,k)
X
≤2 Hm (k(k − 1)P{b1 = m, b2 = m} − (E[ρm (n, k)])2 )1/2
m=3
 n−k+1
X 2
+2 Hm (E[ρm (n, k)])1/2 .
m=3
P 2
n−k+1
As in part (c) we infer m=3 Hm (E[ρm (n, k)])
1/2 = O(n−1 (n − k)2 ) = o(n − k) as n, k → ∞.
Further, for m ≥ 3
(n,k) (n,k)
k(k − 1)P[b1 = m, b2 = m] − (E[ρm (n, k)])2
m−2
k(k + r + s − 1) Y n−k−j
=
n+r+s−1 n+r+s−2−j
j=0
 
m−2 m−2
(k − 1)(k + r + s − 2) Y n−k−m+1−j k(k + r + s − 1) Y n−k−j
× − 
n + r + s − 2m n+r+s−m−1−j n+r+s−1 n+r+s−2−j
j=0 j=0

k(k + r + s − 1)  n−k m−1  n − k − m + 2 m−1

≤
n+r+s−1 n+r+s−2 n+r+s−m
 (k − 1)(k + r + s − 2) k(k + r + s − 1) 
× −
n + r + s − 2m n+r+s−1
2 2
 
k (k + r + s − 1)  n−k  2m−2 1 1
≤ −
n+r+s−1 n+r+s−2 n + r + s − 2m n + r + s − 1
(2m − 1)k 2 (k + r + s − 1)2  n−k 2m−2
=
n + r + s − 2m (n + r + s − 1)2 n + r + s − 2
(2m − 1)k 2  n−k 2m−2
≤ .
n + r + s − 2m n + r + s − 2
Observe that n/(n − 2m) → 1 as n, k → ∞ uniformly in m ∈ N, m ≤ n − k + 1. This in combination
with Lebesgue’s dominated convergence theorem enables us to conclude that

 n−k+1
X  2m − 1 1/2  n−k m−1 2
2
k Hm
n + r + s − 2m n+r+s−2
m=3
(n − k)4  (n − k)2 
∼ 5H32 = o , n, k → ∞.
n3 n
□

Proof of Theorem 7.4. We first show that

hP i
Pn−k+1 Pk (n,j) (n,k) n−k+1 Pk (n,j) (n,k)
m=1 j=1 Z r,s 1[bj = m] − E m=1 j=1 Z r,s 1[bj = m] d
√ −→ N [0; 1] . (7.10)
2 −1 n−k n,k→∞
48
By Lemma 7.5(a) and Chebyshev’s inequality,
hP i
Pk (n,j) (n,k) k (n,j) (n,k)
j=1 Z r,s 1[bj = 1] − E j=1 Z r,s 1[bj = 1] ρ1 (n, k) − E[ρ1 (n, k)] P
√ = √ → 0, n, k → ∞.
n−k n−k
The m = 2 term of the sum gives a principal contribution. Indeed, write
hP i
Pk (n,j) (n,k) k (n,j) (n,k)
j=1 Zr,s 1[bj = 2] − E j=1 Z r,s 1[bj = 2]
√
2 −1 n−k
Pk (n,j) (n,k) r
j=1 (2Zr,s − 3)1[bj = 2] ρ2 (n, k) ρ2 (n, k) − E[ρ2 (n, k)]
= p +3 √ .
ρ2 (n, k) n−k n−k

By Lemma 7.5(b) and Chebyshev’s inequality, the second term converges to 0 in probability. Con-
(n,j) (n,k)
ditionally on Gn,k , kj=1 (2Zr,s − 3)1[bj
P
= 2] has the same distribution as the sum of ρ2 (n, k)
independent random variables taking values ±1 with probability 1/2. Hence, by the central limit
theorem for random walks the first factor converges in distribution to N [0; 1] conditionally on Gn,k ,
hence also unconditionally. By Lemma 7.5 (b),

ρ2 (n, k) P
→ 1, n, k → ∞.
n−k
This secures the distributional convergence of the first term and the sum of two terms to N [0; 1].
Now we prove that the contributions of the counts ρm (n, k), 3 ≤ m ≤ n − k + 1 are negligible.
Write

n−k+1 k h n−k+1 k i
(n,k) (n,k)
X X X X
(n,j) (n,j)
Zr,s 1[bj = m] − E Zr,s 1[bj = m]
m=3 j=1 m=3 j=1
 n−k+1 k n−k+1   n−k+1 
(n,k)
X X X X
(n,j)
= Zr,s 1[bj = m] − Hm ρm (n, k) + Hm (ρm (n, k) − E[ρm (n, k)])
m=3 j=1 m=3 m=3
=: I(n, k) + J(n, k).

We first calculate
n−k+1
X X m  n−k+1
X
E[(I(n, k))2 |Gn,k ] = ℓ−1 (1 − ℓ−1 ) ρm (n, k) ≤ Hm ρm (n, k).
m=3 ℓ=1 m=3

Using now Lemma 7.5(c) we obtain E[I(n, k)2 ] = o(n − k) as n, k → ∞. Finally, by Lemma 7.5(d)
E[J(n, k)2 ] = o(n − k) as n, k → ∞. The proof of (7.10) is complete.
(n,0)
In view of (7.10) and according to Proposition 6.16, it remains to show that E[Zr,s ] = O(1).
But this follows immediately from (6.21). □

Theorems 7.1, 7.2, 7.3 and 7.4 prove Conjecture 2.4 in [24] in a more general form. Hsien-Kuei
Hwang has kindly informed us that the asymptotic normality of the distribution appearing in [24]
was proved by him (unpublished notes) using singularity analysis and the saddle point method.
49
Acknowledgement. We are grateful to Hsien-Kuei Hwang for pointing out the article [24]. ZK
has been supported by the German Research Foundation under Germany’s Excellence Strategy
EXC 2044 - 390685587, Mathematics Münster: Dynamics - Geometry - Structure and by the DFG
priority program SPP 2265 Random Geometric Systems. AM was supported by the Alexander von
Humboldt Foundation.

References
[1] M. Aguiar and S. Mahajan. Topics in hyperplane arrangements. American Mathematical
Society, Providence, RI, 2017.
[2] R. Arratia, A. D. Barbour, and S. Tavaré. Logarithmic combinatorial structures: a probabilistic
approach. European Mathematical Society (EMS), Zürich, 2003.
[3] H. Belbachir and A. Belkhir. Cross recurrence relations for r-Lah numbers. Ars Combinatoria,
110:199–203, 2013.
[4] A. Belkhir. The multivariate Lah and Stirling numbers. J. Integer Seq., 23(4):Art. 20.4.5, 13,
2020.
[5] A. Broder. The r-Stirling numbers. Discrete Mathematics, 49(3):241–259, 1984.
[6] L. Carlitz. Weighted Stirling numbers of the first and second kind. I. Fibonacci Quart.,
18(2):147–162, 1980.
[7] L. Carlitz. Weighted Stirling numbers of the first and second kind. II. Fibonacci Quart.,
18(3):242–257, 1980.
[8] L. Carlitz and R. Scoville. Generalized Eulerian numbers: combinatorial applications. J. Reine
Angew. Math., 265:110–137, 1974.
[9] Ch. A. Charalambides. On a generalized Eulerian distribution. Ann. Inst. Statist. Math.,
43(1):197–206, 1991.
[10] Ch. A. Charalambides. Enumerative combinatorics. Chapman & Hall/CRC, 2002.
[11] Ch. A. Charalambides. Combinatorial methods in discrete distributions. Wiley Ser. Probab.
Stat. John Wiley & Sons, 2005.
[12] Ch. A. Charalambides and M. Koutras. On the differences of the generalized factorials at an
arbitrary point and their combinatorial applications. Discrete Math., 47:183–201, 1983.
[13] Ch. A. Charalambides and J. Singh. A review of the Stirling numbers, their generalizations
and statistical applications. Comm. Statist. Theory Methods, 17(8):2533–2595, 1988.
[14] G.-S. Cheon and J.-H. Jung. r-Whitney numbers of Dowling lattices. Discrete Mathematics,
312:2337–2348, 2012.
[15] H. Crane. The ubiquitous Ewens sampling formula. Statist. Sci., 31(1):1–19, 2016.
[16] L. Devroye. Applications of the theory of records in the study of random trees. Acta Inform.,
26(1-2):123–130, 1988.
[17] M. Dondajewski and J. Szymański. On the distribution of vertex-degrees in a strata of a
random recursive tree. Bull. Acad. Polon. Sci. Sér. Sci. Math., 30(5–6):205–209, 1982.
[18] M. Drmota. Random trees: An interplay between combinatorics and probability. Springer-
Verlag, Vienna, 2009.
[19] A. Gnedin and J. Pitman. Exchangeable Gibbs partitions and Stirling triangles. Zap. Nauchn.
Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov. (POMI), 325:83–102, 244–245, 2005.
[20] T. Godland and Z. Kabluchko. Angle sums of Schläfli orthoschemes. Discrete and Computa-
tional Geometry, 68:125–164, 2022.
50
[21] T. Godland and Z. Kabluchko. Positive hulls of random walks and bridges. Stochastic Processes
and their Applications, 147:327–362, 2022.
[22] R. L. Graham, D. E. Knuth, and O. Patashnik. Concrete mathematics. Addison-Wesley
Publishing Company, second edition, 1994.
[23] A. Gut. Probability: a graduate course. Springer, New York, 2005.
[24] R. van der Hofstad, G. Hooghiemstra, and P. Van Mieghem. Size and weight of shortest path
trees with exponential link weights. Combin. Probab. Comput., 15(6):903–926, 2006.
[25] F. Hoppe. Pólya-like urns and the Ewens’ sampling formula. J. Math. Biol., 20(1):91–94, 1984.
[26] F. Hoppe. Size-biased filtering of Poisson-Dirichlet samples with an application to partition
structures in genetics. J. Appl. Probab., 23(4):1008–1012, 1986.
[27] F. Hoppe. The sampling theory of neutral alleles and an urn model in population genetics. J.
Math. Biol., 25(2):123–159, 1987.
[28] T. Huillet. Occupancy problems related to the generalized Stirling numbers. J. Stat. Phys.,
191(1):Paper No. 5, 2024.
[29] T. Huillet and M. Möhle. On Bernoulli trials with unequal harmonic success probabilities.
Metrika, 2023. to appear.
[30] K. G. Janardan. Relationship between Morisita’s model for estimating the environmental
density and the generalized Eulerian numbers. Ann. Inst. Statist. Math., 40(3):439–450, 1988.
[31] K. G. Janardan. Some properties of the generalized Eulerian distribution. J. Statist. Plann.
Inference, 34(2):159–169, 1993.
[32] S. Janson. Euler-Frobenius numbers and rounding. Online J. Anal. Comb., 8:34, 2013.
[33] K. Joag-Dev and F. Proschan. Negative association of random variables, with applications.
Ann. Statist., 11(1):286–295, 1983.
[34] N. L. Johnson and S. Kotz. Urn models and their application: an approach to modern discrete
probability theory. John Wiley & Sons, Inc., New York, 1977.
[35] N. L. Johnson, S. Kotz, and N. Balakrishnan. Discrete multivariate distributions. John Wiley
& Sons, Inc., New York, 1997.
[36] Z. Kabluchko and A. Marynych. Lah distribution: Stirling numbers, records on compositions,
and convex hulls of high-dimensional random walks. Probability Theory and Related Fields,
184:969–1028, 2022.
[37] Z. Kabluchko, A. Marynych, and H. Pitters. Mod-φ convergence of Stirling distributions and
limit theorems for zeros of their generating functions. J. Math. Anal. Appl., 529(1):Paper No.
127571, 27, 2024.
[38] Z. Kabluchko and D. A. Steigenberger. r-Lah distribution: properties, limit theorems and an
application to compressed sensing. Adv. in Appl. Math., 150:Paper No. 102575, 23, 2023.
[39] Z. Kabluchko, V. Vysotsky, and D. Zaporozhets. Convex hulls of random walks: expected
number of faces and face probabilities. Adv. Math., 320:595–629, 2017.
[40] Z. Kabluchko, V. Vysotsky, and D. Zaporozhets. Convex hulls of random walks, hyperplane
arrangements, and Weyl chambers. Geometric and Functional Analysis, 27(4):880–918, 2017.
[41] J. F. C. Kingman. Poisson processes. The Clarendon Press, Oxford University Press, New
York, 1993.
[42] M. Koutras. Noncentral Stirling numbers and some applications. Discrete Math., 42(1):73–89,
1982.

51
[43] K. Leckey and R. Neininger. Asymptotic analysis of Hoppe trees. J. Appl. Probab., 50(1):228–
238, 2013.
[44] H. M. Mahmoud and R. T. Smythe. On the distribution of leaves in rooted subtrees of recursive
trees. Ann. Appl. Probab., 1(3):406–418, 1991.
[45] R. S. Maier. Triangular recurrences, generalized Eulerian numbers, and related number trian-
gles. Adv. in Appl. Math., 146:Paper No. 102485, 62, 2023.
[46] T. Mansour. Combinatorics of set partitions. CRC Press, Boca Raton, FL, 2013.
[47] T. Mansour and M. Schork. Commutation relations, normal ordering, and Stirling numbers.
CRC Press, Boca Raton, FL, 2016.
[48] I. Mező. The r-Bell numbers. J. Integer Seq., 14(1):Article 11.1.1, 14, 2011.
[49] I. Mező. Recent developments in the theory of Stirling numbers. RIMS Kôkyûroku, 2013:68–80,
2013.
[50] I. Mező. Combinatorics and number theory of counting sequences. CRC Press, 2020.
[51] M. Morisita. Measuring of habitat value by environmental density method. In G. P. et al. Patil,
editor, Statistical Ecology, pages 379–401. The Pennsylvania State University Press, 1971.
[52] D. Najock and C. C. Heyde. On the number of terminal vertices in certain random trees with
an application to stemma construction in philology. J. Appl. Probab., 19(3):675–680, 1982.
[53] K. Nishimura and M. Sibuya. Extended Stirling family of discrete probability distributions.
Comm. Statist. Theory Methods, 26(7):1727–1744, 1997.
[54] G. Nyul and G. Rácz. The r-Lah numbers. Discrete Mathematics, 338(10):1660–1666, 2015.
[55] T. Petersen. Eulerian numbers. Birkhäuser/Springer, New York, 2015.
[56] R. Pinsky. A view from the bridge spanning combinatorics and probability. Enumer. Comb.
Appl., 1(3):Paper No. S2S3, 31, 2021.
[57] J. Pitman. Combinatorial stochastic processes. Springer-Verlag, Berlin, 2006.
[58] R. Shanmugam. On central versus factorial moments. South African Statist. J., 18(2):97–110,
1984.
[59] M. Shattuck. Generalized r-Lah numbers. Proc. Indian Acad. Sci. Math. Sci., 126(4):461–478,
2016.
[60] M. Sibuya. Stirling family of distributions. John Wiley & Sons, 2014.
[61] M. Sibuya and K. Nishimura. Prediction of record-breakings. Statist. Sinica, 7(4):893–906,
1997.
[62] R. T. Smythe and H. M. Mahmoud. A survey of recursive trees. Teor. Ĭmovı̄r. Mat. Stat.,
51:1–29, 1994.
[63] A. J. Stam. Generation of a random partition of a finite set by an urn model. J. Combin.
Theory Ser. A, 35(2):231–240, 1983.
[64] R. Stanley. Enumerative combinatorics. Vol. 1. Cambridge University Press, 1997.
[65] J. Szymański. On the maximum degree and the height of a random recursive tree. In Random
graphs ’87 (Poznań, 1987), pages 313–324. Wiley, Chichester, 1990.
[66] R. van der Hofstad, G. Hooghiemstra, and P. Van Mieghem. First-passage percolation on the
random graph. Probab. Engrg. Inform. Sci., 15(2):225–237, 2001.
[67] P. Van Mieghem. Performance analysis of communications networks and systems. Cambridge
University Press, 2006.
[68] P. Van Mieghem, G. Hooghiemstra, and R. van der Hofstad. A scaling law for the hopcount
in internet. Delft University of Technology, report, 2000125, 2000.

52
[69] P. Van Mieghem, G. Hooghiemstra, and R. van der Hofstad. Stochastic model for the number
of traversed routers in internet. Proceedings of Passive and Active Measurement (PAM), RIPE
NCC (Amsterdam, The Netherlands, 2001), 2001.

53