Mode and Edgeworth Expansion for the Ewens Distribution and the Stirling Numbers

Zakhar Kabluchko Institut für Mathematische Statistik

Universität Münster Orléans–Ring 10 48149 Münster

Germany

zakhar.kabluchko@uni-muenster.de

Alexander Marynych Faculty of Cybernetics

Taras Shevchenko National University of Kyiv 01601 Kyiv Ukraine

marynych@unicyb.kiev.ua

Henning Sulzbach School of Computer Science

McGill University 3480 University Street Montréal, QC H3A 0E9

Canada

henning.sulzbach@gmail.com

Abstract

We provide asymptotic expansions for the Stirling numbers of the first kind and, more generally, the Ewens (or Karamata-Stirling) distribution. Based on these expansions, we obtain some new results on the asymptotic properties of the mode and the maximum of the Stirling numbers and the Ewens distribution. For arbitrary θ > 0 and for all sufficiently large n ∈ ℕ, the unique maximum of the Ewens probability mass function

𝕃_{n}(k) = \frac{θ^{k}}{θ(θ + 1) ⋯ (θ + n − 1)}\left[{n \atop k}\right], k = 1, …, n,

is attained at k = ⌊a_{n}⌋ or ⌈a_{n}⌉, where a_{n} = θ log n − θΓ′(θ)/Γ(θ) − 1/2. We prove that the mode is the nearest integer to a_{n} for a set of n’s of asymptotic density 1, yet this formula is not true for infinitely many n’s.

1 Introduction and statement of results

1.1 Introduction

The (unsigned) Stirling numbers of the first kind \left[{n \atop k}\right] are defined, for n ∈ ℕ and 1 ≤ k ≤ n, by the formula

x^{(n)} := x(x + 1) ⋯ (x + n − 1) = ∑_{k=1}^{n} \left[{n \atop k}\right]x^{k}, x ∈ ℝ.  (1)

For n ∈ ℕ, a random variable K_{n}(θ) is said to have the Ewens distribution with parameter θ > 0 if its probability mass function is given by the formula

ℙ(K_{n}(θ) = k) = \frac{θ^{k}}{θ^{(n)}} \left[{n \atop k}\right], k = 1, …, n.

Bingham [2] called this distribution the Karamata-Stirling law. One can interpret K_{n}(θ) as the number of blocks in a random partition of {1, …, n} distributed according to the Ewens sampling formula, or, equivalently, the number of different alleles in the infinite alleles model, the number of tables in a Chinese restaurant process, or the number of colors in the Hoppe urn. The Ewens sampling formula plays an important role in population genetics [6], [4, Section 1.3]. There is a distributional representation of K_{n}(θ) as a sum of independent random variables

K_{n}(θ) =^{d} ξ_{1} + ⋯ + ξ_{n}, where ξ_{i} ∼ Bern(θ/(θ + i − 1))

and Bern(p) denotes the Bernoulli distribution with parameter p. In the special case θ = 1, classical results going back at least to Feller [7] and Rényi [21] state that the random variable K_{n}(1) has the same distribution as the number of cycles in a uniformly chosen random permutation of n objects, or the number of records in a sample of n i.i.d. variables from a continuous distribution. It follows easily from Lindeberg’s theorem that K_{n}(θ) satisfies a central limit theorem of the form

\frac{K_{n}(θ) − θ log n}{\sqrt{θ log n}} ⟶_{n→∞}^{d} N(0, 1)

known as Goncharov’s CLT in the case θ = 1.

Asymptotic expansions, as n → ∞, of the Stirling numbers \left[{n \atop k}\right] in various regions of k were provided in numerous works [12, 13, 18, 20, 23, 24]. Most notably, Hwang [13, Theorem 2] (and Theorem 14 on page 108 of his dissertation [12] for a more general result) gave an asymptotic expansion valid uniformly in the domain 2 ≤ k ≤ η log n, for any fixed η > 0. Louchard [18, Theorem 2.1] computed three non-trivial terms of the asymptotic expansion in the central regime k = log n + O(\sqrt{log n}) which is similar to the classical Edgeworth expansion in the central limit theorem.

In this short note we start by deriving a full Edgeworth expansion, as n → ∞, for the sequence of probability mass functions k ↦ ℙ(K_{n}(θ) = k) which is uniform both in θ ∈ [1/η, η] (where η > 1) and in k ∈ {1, …, n}; see Theorem 1. Our result is an application of the general Edgeworth expansion for deterministic or random profiles which the authors [16] recently obtained. Using this asymptotic expansion we derive some new results on the mode and the maximum of the Ewens distribution. In the case θ = 1 the mode can be interpreted as the most probable number of cycles in a random permutation of n objects. It was investigated in the works of Hammersley [10, 11] and Erdős [5]. Our results on the mode and the maximum will be stated in Theorems 5 and 7 below.

1.2 Asymptotic expansion of the Ewens distribution

Before stating our main result we need to recall some notions. The (complete) Bell polynomials B_{j}(z_{1}, …, z_{j}) are defined by the formal identity

exp (∑_{j=1}^{∞} \frac{x^{j}}{j!} z_{j}) = ∑_{j=0}^{∞} \frac{x^{j}}{j!} B_{j}(z_{1}, …, z_{j}).

Therefore B_{0} = 1 and, for j ∈ ℕ,

B_{j}(z_{1}, …, z_{j}) = ∑′ \frac{j!}{i_{1}! ⋯ i_{j}!} (\frac{z_{1}}{1!} )^{i_{1}} ⋯ (\frac{z_{j}}{j!} )^{i_{j}} ,  (2)

where the sum ∑′ is taken over all i_{1}, …, i_{j} ∈ ℕ_{0} satisfying 1i_{1} + 2i_{2} + ⋯ + ji_{j} = j. For example, the first three Bell polynomials are given by

B_{1}(z_{1}) = z_{1}, B_{2}(z_{1}, z_{2}) = z_{1}^{2} + z_{2}, B_{3}(z_{1}, z_{2}, z_{3}) = z_{1}^{3} + 3z_{1}z_{2} + z_{3}.  (3)

Further, we will use the “probabilist” Hermite polynomials He_{n}(x) defined by

He_{n}(x) = e^{\frac{1}{2} x^{2}} (−\frac{d}{dx})^{n} e^{− \frac{1}{2} x^{2}}, n ∈ ℕ_{0}.  (4)

The first few Hermite polynomials needed for the first three terms of the expansion are

He_{0}(x) = 1, He_{1}(x) = x, He_{2}(x) = x^{2} − 1, He_{3}(x) = x^{3} − 3x,
He_{4}(x) = x^{4} − 6x^{2} + 3, He_{6}(x) = x^{6} − 15x^{4} + 45x^{2} − 15.

Theorem 1. Fix r ∈ ℕ_{0} and a compact subset L ⊂ (0, ∞). Uniformly over θ ∈ L we have

lim_{n→∞}(log n)^{\frac{r+1}{2}} sup_{k=1,...,n} |ℙ(K_{n}(θ) = k) − \frac{e^{− \frac{1}{2} x_{n}^{2}(k,θ)}}{\sqrt{2πθ log n}} ∑_{j=0}^{r} \frac{H_{j}(x_{n}(k, θ))}{(θ log n)^{j/2}} | = 0.

Here, x_{n}(k, θ) = \frac{k−θ log n}{\sqrt{θ log n}} and H_{j}(x) is a polynomial of degree 3j given by

H_{j}(x) := H_{j}(x, θ) = \frac{(−1)^{j}}{j!} e^{\frac{1}{2} x^{2}}B_{j}(D̃_{1}, …, D̃_{j})e^{− \frac{1}{2} x^{2}},  (5)

where B_{j} is the j-th Bell polynomial and D̃_{1}, D̃_{2}, … are differential operators given by

D̃_{j} := D̃_{j}(θ) = \frac{1}{(j + 1)(j + 2)} (\frac{d}{dx})^{j+2} + χ̃_{j}(0) (\frac{d}{dx})^{j}  (6)

with χ̃_{j}(β) = − (\frac{d}{dβ} )^{j} log Γ(θe^{β}) and Γ denoting the Euler gamma function.

Remark 2. It follows from (3), (5) and (6) that the first three coefficients of the expansion are

H_{0}(x) = 1,
H_{1}(x) = −\frac{Γ′(θ)}{Γ(θ)}θx + \frac{1}{6}He_{3}(x),
H_{2}(x) = (θ^{2} \frac{Γ^{′2}(θ)}{Γ^{2}(θ)} − \frac{θ^{2}Γ′′(θ) + θΓ′(θ)}{2Γ(θ)}) He_{2}(x) + (\frac{1}{24} − \frac{Γ′(θ)}{6Γ(θ)}θ) He_{4}(x) + \frac{1}{72}He_{6}(x).

An expression for χ̃_{j}(0) involving polygamma functions and Stirling numbers of the second kind will be given below in (13). The tilde in D̃_{j} and χ̃_{j} is needed to keep the notation consistent with our more general work [16]. It is easy to check that H_{j}(−x) = (−1)^{j}H_{j}(x) [16, Remark 2.4].

To compute H_{j}(x) one can proceed as follows. First, express \frac{1}{j!}B_{j}(D̃_{1}, …, D̃_{j}) as a polynomial in D := \frac{d}{dx} (and note that only even/odd powers of D are present if j is even/odd). Then replace each occurrence of D^{l} by He_{l}(x); see (4) for justification.

Remark 3. It is possible to choose the value of θ as a function of k. One natural choice is θ = 1 which provides a full version of Louchard’s expansion [18, Theorem 2.1] (although he used a slightly different normalization in his analogue of x_{n}(k, 1) and his term −355x^{3}/144 should be replaced by −47x^{3}/144). Another possible choice is θ = k/ log n (so that x_{n}(k, θ) = 0), which gives a large-deviation-type expansion valid uniformly in the region η^{−1} log n < k < η log n, for fixed η > 1 and q ∈ ℕ_{0}:

\frac{(k/ log n)^{k}}{(k/ log n)^{(n)}} \left[{n \atop k}\right] = \frac{1}{\sqrt{2πk}} ∑_{s=0}^{q} \frac{H_{2s}(0, k/ log n)}{k^{s}} + o (\frac{1}{(log n)^{q+1}} ) .

Observe that the terms with half-integer powers of k are not present in the sum because H_{2j+1}(0) = 0. Using the formula

\frac{Γ(n + θ)}{n!} = n^{θ−1} (1 + O (\frac{1}{n}))

yields the expansion

\frac{1}{n!}\left[{n \atop k}\right] = \frac{1}{Γ(θ)}n^{θ−θ log θ−1} (\frac{1}{\sqrt{2πk}} ∑_{s=0}^{q} \frac{H_{2s}(0, θ)}{k^{s}} + o (\frac{1}{(log n)^{q+1}} ))  (7)

valid as n → ∞ uniformly over k in the region θ = k/ log n ∈ (η^{−1}, η). In this region, this expansion must be equivalent to Hwang’s result [13, Theorem 2]. It is not easy to rigorously verify this equivalence by a direct comparison, but we checked using Mathematica 9 that the first three non-trivial terms coincide. Note a misprint in the formula for the remainder term Z_{μ}(m, n) in Hwang [13, Theorem 2]: (log n)^{m}/(m!n) should be replaced by (log n)/(mn). Expansion (7) could be also deduced from the work of Féray et al. [8, Theorem 3.4].

Taking sums over k in Theorem 1 and using the Euler-Maclaurin formula to approximate Riemann sums by integrals, one obtains that

ℙ (\frac{K_{n}(θ) − θ log n}{\sqrt{θ log n}} ≤ x) = Φ(x) + \frac{e^{−x^{2}/2}}{\sqrt{2πθ log n}} (\frac{1}{2} − \frac{x^{2} − 1}{6} + θ\frac{Γ′(θ)}{Γ(θ)}) + O (\frac{1}{log n}) ,

uniformly in x ∈ (θ log n)^{−1/2}(ℤ − θ log n), where Φ(x) is the standard normal distribution function. The proof follows Grübel and Kabluchko [9, Proposition 2.5] and is therefore omitted. Yamato [25] recently stated a slightly incorrect version of this expansion missing the term 1/2 which comes from the Euler-Maclaurin formula. Similarly, one can obtain further terms in the expansion of the distribution function of (K_{n}(θ) − θ log n)/\sqrt{θ log n}.

Remark 4. Since the set L ⊆ (0, ∞) in Theorem 1 has to be chosen compact, our results do not yield asymptotic expansions for ℙ(K_{n}(θ) = k) in the regime k = o(log n) of the same precision as Hwang’s [13, Theorems 1 and 2]. Also, they do not extend straightforwardly to the case k = n − O(n^{α}) for 0 < α < 1 treated by Louchard [18, Section 3]. A generalization of our approach to cover these regions will be content of future work.

1.3 Mode and maximum of the Ewens distribution

Theorem 1 allows us to deduce various results on the mode and the maximum of the Ewens distribution. A mode is any value k ∈ {1, …, n} maximizing ℙ(K_{n}(θ) = k), while the maximum M_{n}(θ) is defined by

M_{n}(θ) = max_{1≤k≤n} ℙ(K_{n}(θ) = k).

Let u_{n}(θ) denote the least mode. In this context, it is important to note that, for all θ > 0, the function k ↦ ℙ(K_{n}(θ) = k) is log-concave by a theorem attributed to Newton [11, 22], and

ℙ(K_{n}(θ) = 1) < … < ℙ(K_{n}(θ) = u_{n}(θ))
≥ ℙ(K_{n}(θ) = u_{n}(θ) + 1) > … > ℙ(K_{n}(θ) = n).  (8)

In particular, there are at most two modes. For θ = 1, Erdős [5], proving a conjecture of Hammersley [11], showed that the mode is unique for all n ≥ 3. By (8), uniqueness also holds for irrational θ; however, for rational θ, the mode need not be unique since, for example,

\frac{2}{3}\left[{3 \atop 1}\right] = (\frac{2}{3})^{2} \left[{3 \atop 2}\right] > (\frac{2}{3})^{3} \left[{3 \atop 3}\right].

Theorem 5. Fix θ > 0. There exists N_{1} ∈ ℕ such that for n ≥ N_{1}, u_{n}(θ) is the unique mode of the Ewens distribution with parameter θ. The mode u_{n}(θ) equals one of the numbers ⌊u_{n}^{∗}(θ)⌋ or ⌈u_{n}^{∗}(θ)⌉, where

u_{n}^{∗}(θ) = θ log n − \frac{θΓ′(θ)}{Γ(θ)} − \frac{1}{2}

and ⌊·⌋, ⌈·⌉ denote the floor and the ceiling functions, respectively. Write δ_{n}(θ) := min_{k∈ℤ} |u_{n}^{∗}(θ) − k|. For the maximum M_{n}(θ), we have

\sqrt{2πθ log n} M_{n}(θ) = 1 + \frac{θ(log Γ)′(θ) + θ^{2}(log Γ)′′(θ) + 1/12 − δ_{n}^{2}(θ)}{2θ log n} + o (\frac{1}{log n}) .

In the case θ = 1, Hammersley [11] and Erdős [5] derived related results for the mode. Cramer [3] discusses statistical applications and Mező [19] provides an overview and generalizations. Theorem 5 states that the mode is one of the numbers ⌊log n + γ − \frac{1}{2}⌋ or ⌈log n + γ − \frac{1}{2}⌉, for sufficiently large n. In fact, this holds for all n ∈ ℕ.

Proposition 6. u_{n}(1) ∈ \{⌊log n + γ − \frac{1}{2}⌋, ⌈log n + γ − \frac{1}{2}⌉\} for all n ∈ ℕ.

The proof uses the following formula of Hammersley [11]:

u_{n}(1) = ⌊log n + γ + \frac{ζ(2) − ζ(3)}{log n + γ − \frac{3}{2}} + \frac{h(n)}{(log n + γ − \frac{3}{2})^{2}} ⌋ ,  (9)

for some −1.098011 < h(n) < 1.430089. Hwang [12, Section 5.7.9] gives a more precise expansion. Erdős [5] observed that, for n > 189, Hammersley’s formula implies that the mode is one of the numbers ⌊log(n − 1) + \frac{1}{2}⌋ or ⌊log(n − 1) + 1⌋. Note that his Σ_{n,s} equals \left[{n+1 \atop n+1−s}\right] and his n − f(n) is u_{n+1}(1) − 1 in our notation.

The next theorem provides more precise information about the behavior of the mode. Recall that a set A ⊂ ℕ has asymptotic density α ∈ [0, 1] if

lim_{n→∞} \frac{#(A ∩ \{1, …, n\})}{n} = α.

For x ∈ ℝ, let {x} = x − ⌊x⌋ denote the fractional part of x. Let nint(x) be the integer closest to x (if {x} = 1/2, we agree to take nint(x) = ⌈x⌉). That is,

nint(x) := arg min_{k∈ℤ} |x − k| = ⌊x + \frac{1}{2}⌋ .

Theorem 7. Fix θ > 0. The mode u_{n}(θ) of the Ewens distribution with parameter θ has the following properties:

(i) there exists a constant C_{0} > 0 such that, for all n ∈ ℕ satisfying

|\{u_{n}^{∗}(θ)\} − \frac{1}{2}| > \frac{C_{0}}{log n},

the mode u_{n}(θ) equals

nint(u_{n}^{∗}(θ)) = ⌊θ log n − \frac{θΓ′(θ)}{Γ(θ)}⌋;

(ii) there are arbitrarily long intervals of consecutive n’s for which u_{n}(θ) = ⌈u_{n}^{∗}(θ)⌉; similarly, there are arbitrarily long intervals of consecutive n’s for which u_{n}(θ) = ⌊u_{n}^{∗}(θ)⌋;

(iii) the set of n ∈ ℕ such that u_{n}(θ) = nint(u_{n}^{∗}(θ)) has asymptotic density one;

(iv) there are infinitely many n ∈ ℕ such that u_{n}(θ) ≠ nint(u_{n}^{∗}(θ)).

The proof of part (iv) uses five terms in the Edgeworth expansion, where the first two terms influence the form of u_{n}^{∗}(θ), while the remaining terms are needed for technical reasons. The idea is that the formula u_{n}(θ) = nint(u_{n}^{∗}(θ)) becomes wrong if the fractional part of u_{n}^{∗}(θ) is slightly below \frac{1}{2}, so that higher order terms in the Edgeworth expansion decide which of the values ⌊u_{n}^{∗}(θ)⌋ and ⌈u_{n}^{∗}(θ)⌉ is the mode. Using even more terms in the expansion, it is possible to replace u_{n}^{∗}(θ) by some more complicated expressions involving higher-order corrections in inverse powers of θ log n [12, Section 5.7.9]. However, it seems that there is no formula of the form

u_{n}(1) = nint (log n + a_{0} + \frac{a_{1}}{log n} + ⋯ + \frac{a_{r}}{(log n)^{r}} )

which is valid for all sufficiently large n.

Finally, we would like to mention that one can easily obtain counterparts of the above results for the B- and D-analogues of Stirling numbers of the first kind. These are defined as the coefficients of (x + 1)(x + 3) ⋯ (x + 2n − 1) and ((x + 1)(x + 3) ⋯ (x + 2n − 3))(x + n − 1), respectively. They appear, for example, in the study of intrinsic volumes of Weyl chambers [14].

2 Proofs

Proof of Theorem 1. The proof follows from the general Edgeworth expansion for random or deterministic profiles [16, Theorem 2.1]. We consider the sequence of “profiles”

𝕃_{n}(k) := ℙ(K_{n}(θ) = k) = \frac{θ^{k}}{θ^{(n)}} \left[{n \atop k}\right]𝟙_{\{k∈\{1,...,n\}\}},

and define

w_{n} := θ log n, φ(β) := e^{β} − 1, (β_{−}, β_{+}) = ℝ, 𝒟 = {z ∈ ℂ: |Im z| < π}.

In order to apply [16, Theorem 2.1], we need to check Conditions A1–A4 given in the beginning of Section 2 of the cited paper. Note that

W_{n}(β) := e^{−φ(β)w_{n}} ∑_{k∈ℤ} e^{βk}𝕃_{n}(k) = n^{−θ(e^{β}−1)} ∑_{k=1}^{n} e^{βk} \frac{θ^{k}}{θ^{(n)}} \left[{n \atop k}\right]

= n^{−θ(e^{β}−1)} \frac{(θe^{β})^{(n)}}{θ^{(n)}} = n^{−θ(e^{β}−1)} \frac{Γ(θe^{β} + n)Γ(θ)}{Γ(θe^{β})Γ(θ + n)} ⟶_{n→∞} \frac{Γ(θ)}{Γ(θe^{β})} =: W_{∞}(β)

locally uniformly in β ∈ 𝒟 with a rate of convergence which is polynomial in n^{−1}. Hence Conditions A1–A3 are satisfied. In order to check A4 it is enough to show that for every a > 0, r ∈ ℕ and every compact subset K_{1} of ℝ

sup_{β∈K_{1}} sup_{a≤u≤π} (n^{−θ(e^{β}−1)} |\frac{Γ(θe^{β+iu} + n)Γ(θ)}{Γ(θ + n)Γ(θe^{β+iu})}|) = o(log^{−r} n), n → ∞.

But this easily follows from

sup_{β∈K_{1}} (n^{−θ(e^{β}−1)} sup_{a≤u≤π} |\frac{Γ(θe^{β+iu} + n)Γ(θ)}{Γ(θ + n)Γ(θe^{β+iu})}|)

≤ C sup_{β∈K_{1}} (n^{−θ(e^{β}−1)} sup_{a≤u≤π} |\frac{Γ(θe^{β+iu} + n)}{Γ(θ + n)}|) ≤ C_{1} sup_{β∈K_{1}} n^{θe^{β}(cos a−1)},

with constants C, C_{1} depending on K_{1}, θ and a. Therefore, Theorem 2.1 of [16] is applicable for the Ewens distribution with arbitrary fixed θ > 0. In particular, for θ = 1, we obtain

(log n)^{\frac{r+1}{2}} sup_{β∈K} sup_{1≤k≤n} |\frac{Γ(e^{β})e^{βk}}{n^{e^{β}−1}n!} \left[{n \atop k}\right] − \frac{e^{− \frac{1}{2} x_{n}^{2}(k,e^{β})}}{\sqrt{2πe^{β} log n}} ∑_{j=0}^{r} \frac{G_{j}(x_{n}(k, e^{β}); β)}{(log n)^{j/2}} | ⟶_{n→∞} 0,  (10)

where K is a compact subset of ℝ, and the polynomials G_{0}, G_{1}, … are defined as in Theorem 2.1 of [16]: for j ∈ ℕ_{0}, we have

G_{j}(x; β) = \frac{(−1)^{j}}{j!} e^{\frac{1}{2} x^{2}}B_{j}(D_{1}, …, D_{j})e^{− \frac{1}{2} x^{2}}  (11)

with the differential operators

D_{j} := D_{j}(β) = e^{− \frac{1}{2} βj} (\frac{1}{(j + 1)(j + 2)} (\frac{d}{dx})^{j+2} + χ_{j}(β) (\frac{d}{dx})^{j}) ,  (12)

where

χ_{j}(β) = − (\frac{d}{dβ} )^{j} log Γ(e^{β}).

Now, if L ⊆ (0, ∞) is compact, then K := log L is compact in ℝ. Applying (10) with K = log L and β = log θ ∈ K, we obtain

(log n)^{\frac{r+1}{2}} sup_{θ∈L} sup_{1≤k≤n} |\frac{Γ(θ)θ^{k}}{n^{θ−1}n!}\left[{n \atop k}\right] − \frac{e^{− \frac{1}{2} x_{n}^{2}(k,θ)}}{\sqrt{2πθ log n}} ∑_{j=0}^{r} \frac{G_{j}(x_{n}(k, θ); log θ)}{(log n)^{j/2}} | ⟶_{n→∞} 0.

By Stirling’s formula, uniformly in θ ∈ L, n ∈ ℕ and 1 ≤ k ≤ n, we have

\frac{Γ(θ)θ^{k}}{n^{θ−1}n!}\left[{n \atop k}\right] = \frac{θ^{k}}{θ^{(n)}} \left[{n \atop k}\right](1 + O(n^{−1})) = \frac{θ^{k}}{θ^{(n)}} \left[{n \atop k}\right] + O(n^{−1}).

We conclude the proof by noting that G_{j}(x; log θ) = θ^{−j/2}H_{j}(x) which follows directly from χ̃_{j}(0) = χ_{j}(log θ). Indeed, by comparing (6) and (12), we obtain

D_{j}(log θ) = θ^{−j/2}D̃_{j}(θ),

which implies that

B_{j}(D_{1}(log θ), …, D_{j}(log θ)) = θ^{−j/2}B_{j}(D̃_{1}(θ), …, D̃_{j}(θ))

since B_{j}(z_{1}, …, z_{j}) is a sum of terms of the form c · z_{1}^{i_{1}}z_{2}^{i_{2}} ⋯ z_{j}^{i_{j}} with 1i_{1} + 2i_{2} + ⋯ + ji_{j} = j; see (2). Comparing (5) and (11), we obtain the required identity G_{j}(x; log θ) = θ^{−j/2}H_{j}(x).

To see that χ̃_{j}(0) = χ_{j}(log θ), one can easily show by induction over j ≥ 1 that, both

χ_{j}(β) = − ∑_{ℓ=1}^{j} \left\{{j \atop ℓ}\right\}ψ^{(ℓ−1)}(e^{β})e^{ℓβ},

and

χ̃_{j}(β) = − ∑_{ℓ=1}^{j} \left\{{j \atop ℓ}\right\}ψ^{(ℓ−1)}(θe^{β})(θe^{β})^{ℓ}.  (13)

Here ψ^{(j)}(x) = (log Γ(x))^{(j+1)} denotes the polygamma function and \left\{{n \atop k}\right\} is the Stirling number of the second kind satisfying the recurrence

\left\{{n + 1 \atop k}\right\} = \left\{{n \atop k − 1}\right\} + k\left\{{n \atop k}\right\}, 1 ≤ k ≤ n, n ∈ ℕ,

with initial conditions \left\{{0 \atop 0}\right\} = 1, \left\{{n \atop 0}\right\} = \left\{{0 \atop n}\right\} = 0.

Proof of Theorem 5. It follows from Theorems 2.10 in [16] that for sufficiently large n, the maximizers of the function k ↦ ℙ(K_{n}(θ) = k) must be of the form ⌊u_{n}^{∗}⌋ or ⌈u_{n}^{∗}⌉.

Next we prove that the maximizer is unique (for sufficiently large n) by following a method of Erdős [5] who considered the case θ = 1. Thanks to (8), the uniqueness is evident if θ is irrational. Hence, we assume that θ = Q_{1}/Q_{2} is rational with Q_{1}, Q_{2} being integer. We have, by (1),

\left[{n \atop k}\right] = ∑_{1≤a_{1}<...<a_{n−k}≤n−1} a_{1} ⋯ a_{n−k}.

Put k_{n} = ⌈u_{n}^{∗}(θ)⌉ = θ log n + O(1) as n → ∞. By (8), it is sufficient to show that

θ^{k_{n}} \left[{n \atop k_{n}}\right] ≠ θ^{k_{n}−1}\left[{n \atop k_{n} − 1}\right].  (14)

By Erdős’ argument relying on the prime number theorem with an appropriate error term [5, p. 233], for all sufficiently large n, there is a prime number p satisfying (n − 1)/k_{n} < p < (n − 1)/(k_{n} − 1). Then,

\left[{n \atop k_{n}}\right] ≢ 0 (mod p), \left[{n \atop k_{n} − 1}\right] ≡ 0 (mod p)

because in the representation of the former Stirling number all products except one are divisible by p, whereas in the latter all products are divisible by p. If n is large, p is not among the prime factors of Q_{1} and Q_{2}. Hence (14) follows and the mode of K_{n}(θ) is unique. Finally, the formula for M_{n} follows from Theorem 2.13 of [16].

Proof of Proposition 6. Recall Hammersley’s formula (9):

u_{n}(1) = ⌊log n + γ + \frac{ζ(2) − ζ(3)}{log n + γ − \frac{3}{2}} + \frac{h(n)}{(log n + γ − \frac{3}{2})^{2}} ⌋

with some −1.1 < h(n) < 1.44. It is easy to check that

\frac{ζ(2) − ζ(3)}{x} − \frac{1.1}{x^{2}} > −\frac{1}{2} and \frac{ζ(2) − ζ(3)}{x} + \frac{1.44}{x^{2}} < \frac{1}{2}

for x > 2.5. Hence, the proposition is true for log n + γ − \frac{3}{2} > 2.5, that is for n ≥ 31. For n = 1, 2, …, 30 the statement is easy to verify using Mathematica 9.

Proof of Theorem 7 (i) and (ii). Part (i) follows essentially from Theorem 2.10 in [16] and its proof. Namely, by [16, Equation (90)], for k = k(n) = u_{n}^{∗}(θ) + g ∈ ℤ with g = O(1), we have

\sqrt{2πθ log n} (ℙ(K_{n}(θ) = k + 1) − ℙ(K_{n}(θ) = k)) = −\frac{2g + 1}{2θ log n} + o (\frac{1}{log n}) .

The same relation, but with a better remainder term O(\frac{1}{log^{2} n}), follows from (16) which we shall prove below. Taking g = −{u_{n}^{∗}(θ)}, so that k = ⌊u_{n}^{∗}(θ)⌋ and k + 1 = ⌈u_{n}^{∗}(θ)⌉, yields

ℙ(K_{n}(θ) = ⌈u_{n}^{∗}(θ)⌉) − ℙ(K_{n}(θ) = ⌊u_{n}^{∗}(θ)⌋) = \frac{1}{\sqrt{2πθ log n}} (\frac{\{u_{n}^{∗}(θ)\} − \frac{1}{2}}{θ log n} + O (\frac{1}{log^{2} n})) .

It follows that there is a sufficiently large constant C_{0} > 0 such that, if \{u_{n}^{∗}(θ)\} > \frac{1}{2} + \frac{C_{0}}{log n}, then the right-hand side is positive, and the mode equals ⌈u_{n}^{∗}(θ)⌉. Similarly, if {u_{n}^{∗}(θ)} < \frac{1}{2} − \frac{C_{0}}{log n}, then the right-hand side is negative, and the mode equals ⌊u_{n}^{∗}(θ)⌋.

The proof of part (ii) follows immediately from part (i) and the fact that, for every fixed L > 0, we have log(n + L) − log n → 0 as n → ∞.

Proof of Theorem 7 (iii). In view of part (i) it suffices to show that

lim sup_{ε→0} lim sup_{n→∞} \frac{#\{1 ≤ k ≤ n : dist(u_{k}^{∗}(θ), ℤ + 1/2) < ε\}}{n} = 0,

which, in turn, follows from the fact that

lim sup_{ε→0} lim sup_{n→∞} \frac{#\{1 ≤ k ≤ n : dist(log k, αℤ + β) < ε\}}{n} = 0,  (15)

for all α > 0 and β ∈ ℝ. Equation (15) would be true if the sequence of fractional parts of α^{−1} log k, k ∈ ℕ, were uniformly distributed on [0, 1]. However, the latter claim is unfortunately not true [17, Examples 2.4 and 2.5, pp. 8–9]. Let us prove (15). We have, assuming that ε < α/2,

#{1 ≤ k ≤ n : dist(log k, αℤ + β) < ε} = ∑_{k=1}^{n} #{j ∈ ℤ : dist(log k, αj + β) < ε}

= ∑_{j∈ℤ} #{1 ≤ k ≤ n : e^{αj+β−ε} < k < e^{αj+β+ε}}

≤ ∑_{j∈ℤ} #{k ∈ ℕ : e^{αj+β−ε} ∨ 1 ≤ k ≤ e^{αj+β+ε} ∧ n}.

The summand on the right-hand side is the number of integers in the interval [e^{αj+β−ε} ∨ 1, e^{αj+β+ε} ∧ n] (which is empty if either e^{αj+β−ε} > n or e^{αj+β+ε} < 1). Hence, it is bounded from above by (e^{αj+β+ε} ∧ n − e^{αj+β−ε} ∨ 1 + 1)_{+}. Therefore,

#{1 ≤ k ≤ n : dist(log k, αℤ + β) < ε} ≤ ∑_{j∈ℤ} (e^{αj+β+ε} ∧ n − e^{αj+β−ε} ∨ 1 + 1)_{+} .

Further,

∑_{j∈ℤ} (e^{αj+β+ε} ∧ n − e^{αj+β−ε} ∨ 1 + 1)_{+}

= ∑_{j∈ℤ} e^{αj+β+ε}𝟙_{{αj+β+ε<0}} + ∑_{j∈ℤ} e^{αj+β+ε}𝟙_{{αj+β−ε<0,0≤αj+β+ε<log n}}

+ ∑_{j∈ℤ} n𝟙_{{αj+β−ε<0,log n≤αj+β+ε}}

+ ∑_{j∈ℤ} (e^{αj+β+ε} − e^{αj+β−ε} + 1) 𝟙_{{αj+β−ε≥0,αj+β+ε<log n}}

+ ∑_{j∈ℤ} (n − e^{αj+β−ε} + 1)_{+} 𝟙_{{αj+β−ε≥0,log n≤αj+β+ε}}.

Note that the first series converges, the second contains at most one summand since we assume ε < α/2, and the third vanishes for n large enough. It can be checked that

∑_{j∈ℤ} (e^{αj+β+ε} − e^{αj+β−ε} + 1) 𝟙_{{αj+β−ε≥0,αj+β+ε<log n}} ≤ C(α, β)(e^{β+ε} − e^{β−ε})n

with an absolute constant C(α, β). Further, for n sufficiently large, we have

∑_{j∈ℤ} (n − e^{αj+β−ε} + 1)_{+} 𝟙_{{αj+β−ε≥0,log n≤αj+β+ε}} ≤ n(1 − e^{−2ε}) + 1.

Putting pieces together gives (15).

Proof of Theorem 7 (iv). Recall the notation w_{n} = θ log n and x_{n}(k) = x_{n}(k, θ) = (k − w_{n})/\sqrt{w_{n}}. Using Theorem 1 with r = 4, we obtain

\sqrt{2πw_{n}} ℙ(K_{n}(θ) = k) = e^{− \frac{1}{2} x_{n}^{2}(k)}

× (1 + \frac{H_{1}(x_{n}(k))}{w_{n}^{1/2}} + \frac{H_{2}(x_{n}(k))}{w_{n}} + \frac{H_{3}(x_{n}(k))}{w_{n}^{3/2}} + \frac{H_{4}(x_{n}(k))}{w_{n}^{2}} + o (\frac{1}{log^{2} n})) ,

as n → ∞ uniformly in 1 ≤ k ≤ n. Now let k = θ log n + a, where a = O(1) as n → ∞, so that x_{n}(k) = a/w_{n}^{1/2}. We have

H_{1}(x_{n}(k)) = A_{11}(θ)\frac{a}{w_{n}^{1/2}} + A_{12}(θ)\frac{a^{3}}{w_{n}^{3/2}} ,
H_{2}(x_{n}(k)) = A_{21}(θ) + A_{22}(θ)\frac{a^{2}}{w_{n}} + o (\frac{1}{w_{n}} ) ,
H_{3}(x_{n}(k)) = A_{31}(θ)\frac{a}{w_{n}^{1/2}} + o (\frac{1}{w_{n}^{1/2}} ) ,
H_{4}(x_{n}(k)) = A_{41}(θ) + o(1),

where A_{11}(θ), …, A_{41}(θ) are some polynomials in χ̃_{1}(0), χ̃_{2}(0), χ̃_{3}(0) and χ̃_{4}(0); see Remark 2. Plugging these expressions into the asymptotic expansion above and using the expansion e^{y} = 1 + y + y^{2}/2 + o(y^{2}), as y → 0, yields

\sqrt{2πw_{n}} ℙ(K_{n}(θ) = k) = 1 − (\frac{a^{2}}{2} − A_{11}(θ)a − A_{21}(θ)) \frac{1}{w_{n}} + \frac{P_{θ}(a)}{w_{n}^{2}} + o (\frac{1}{log^{2} n}) ,

where

P_{θ}(a) := \frac{1}{8}a^{4} + (A_{12}(θ) − \frac{1}{2}A_{11}(θ)) a^{3} + (A_{22}(θ) − \frac{1}{2}A_{21}(θ)) a^{2} + A_{31}(θ)a + A_{41}(θ).

Now let us write k = θ log n + a^{∗} + g, where a^{∗} := A_{11}(θ) = −\frac{θΓ′(θ)}{Γ(θ)} − \frac{1}{2}, yielding

\sqrt{2πw_{n}} ℙ(K_{n}(θ) = k) = 1 − (\frac{g^{2} − (a^{∗})^{2}}{2} − A_{21}(θ)) \frac{1}{w_{n}} + \frac{P_{θ}(a^{∗} + g)}{w_{n}^{2}} + o (\frac{1}{log^{2} n}) . (16)

We are interested in g being either ⌊u_{n}^{∗}(θ)⌋ − u_{n}^{∗}(θ) =: g′_{n} or ⌈u_{n}^{∗}(θ)⌉ − u_{n}^{∗}(θ) =: g′′_{n}. Let M be the set of natural numbers n with {u_{n}^{∗}(θ)} < 1/2 < {u_{n+1}^{∗}(θ)}. Note that M has infinitely many elements because log n → ∞ and log(n + 1) − log n → 0. In the remainder of the proof, we always consider n ∈ M. Since u_{n+1}^{∗}(θ) − u_{n}^{∗}(θ) = O(n^{−1}), we have

g′_{n} = −1/2 + O(n^{−1}), g′′_{n} = 1/2 + O(n^{−1}).

Putting k = ⌊u_{n}^{∗}(θ)⌋ into (16) yields

\sqrt{2πw_{n}} ℙ(K_{n}(θ) = ⌊u_{n}^{∗}(θ)⌋)

= 1 − (\frac{(g′_{n})^{2} − (a^{∗})^{2}}{2} − A_{21}(θ)) \frac{1}{w_{n}} + \frac{P_{θ}(a^{∗} + g′_{n})}{w_{n}^{2}} + o (\frac{1}{log^{2} n})

= 1 − (\frac{1 − 4(a^{∗})^{2}}{8} − A_{21}(θ)) \frac{1}{w_{n}} + \frac{P_{θ}(a^{∗} − 1/2)}{w_{n}^{2}} + o (\frac{1}{log^{2} n}) .

Analogously, putting k = ⌈u_{n}^{∗}(θ)⌉ gives

\sqrt{2πw_{n}} ℙ(K_{n}(θ) = ⌈u_{n}^{∗}(θ)⌉) = 1 − (\frac{1 − 4(a^{∗})^{2}}{8} − A_{21}(θ)) \frac{1}{w_{n}} + \frac{P_{θ}(a^{∗} + 1/2)}{w_{n}^{2}} + o (\frac{1}{log^{2} n}) .

For sufficiently large n the mode u_{n}(θ) equals either ⌊u_{n}^{∗}(θ)⌋ or ⌈u_{n}^{∗}(θ)⌉ depending on the sign of

s^{∗}(θ) := P_{θ}(a^{∗} + 1/2) − P_{θ}(a^{∗} − 1/2).

In the following we shall show that s^{∗}(θ) > 0, hence u_{n}(θ) = ⌈u_{n}^{∗}(θ)⌉, while nint(u_{n}(θ)) = ⌊u_{n}^{∗}(θ)⌋, so that u_{n}(θ) ≠ nint(u_{n}^{∗}(θ)). Recalling the polygamma function ψ^{(m)}(θ) = (log Γ(θ))^{(m+1)}, the authors [15] checked with the help of Mathematica 9 that

s^{∗}(θ) = \frac{θ^{2}}{2} (2ψ^{(1)}(θ) + θψ^{(2)}(θ)) .

Using the well-known formula for the polygamma function [1, 6.4.10]

ψ^{(m)}(θ) = (log Γ(θ))^{(m+1)} = (−1)^{m+1}m! ∑_{k=0}^{∞} \frac{1}{(θ + k)^{m+1}} , −θ ∉ ℕ_{0}, m ≥ 1,

we finally obtain

s^{∗}(θ) = θ^{2} ∑_{k=1}^{∞} \frac{k}{(θ + k)^{3}} , θ > 0,

yielding positivity of s^{∗}(θ) for all θ > 0. The proof of part (iv), as well as of the whole theorem, is complete.

Remark 8. For θ = 1 we have s^{∗}(1) = ζ(2) − ζ(3), a term appearing in Hammersley’s formula (9). In fact, in the special case θ = 1 part (iv) could be deduced directly from (9).

3 Acknowledgements

The work of Alexander Marynych was supported by a Humboldt Research Fellowship of the Alexander von Humboldt Foundation. The work of Henning Sulzbach was supported by a Feodor Lynen Research Fellowship of the Alexander von Humboldt Foundation.

References

[1] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, volume 55 of National Bureau of Standards Applied Mathematics Series, 1964.

[2] N. H. Bingham, Tauberian theorems for Jakimovski and Karamata-Stirling methods, Mathematika 35 (1988), 216–224.

[3] E. Cramer, Asymptotic estimators of the sample size in a record model, Statist. Papers 41 (2000), 159–171.

[4] R. Durrett, Probability Models for DNA Sequence Evolution, Probability and its Applications, Springer, second edition, 2008.

[5] P. Erdős, On a conjecture of Hammersley, J. London Math. Soc. 28 (1953), 232–236.

[6] W. Ewens and S. Tavaré, Multivariate Ewens distribution. In N. L. Johnson, S. Kotz, and N. Balakrishnan, editors, Discrete Multivariate Distributions, Wiley Series in Probability and Statistics: Applied Probability and Statistics. John Wiley & Sons, 1997.

[7] W. Feller, The fundamental limit theorems in probability, Bull. Amer. Math. Soc. 51 (1945), 800–832.

[8] V. Féray, P.-L. Méliot, and A. Nikeghbali, Mod–ϕ convergence, I: Normality zones and precise deviations, preprint, 2015, http://arxiv.org/abs/1304.2934.

[9] R. Grübel and Z. Kabluchko, Edgeworth expansions for profiles of lattice branching random walks, to appear in Ann. Inst. H. Poincaré, 2015,

http://arxiv.org/abs/1503.04616.

[10] J. M. Hammersley, On estimating restricted parameters, J. Roy. Statist. Soc. Ser. B. 12 (1950), 192–229, discussion, 230–240.

[11] J. M. Hammersley, The sums of products of the natural numbers, Proc. London Math. Soc. (3) 1 (1951), 435–452.

[12] H.-K. Hwang, Théorèmes limites pour les structures combinatories et les fonctions

arithmétiques, Ph.D. thesis, 1994, http://algo.stat.sinica.edu.tw/hk/?p=80.

[13] H.-K. Hwang, Asymptotic expansions for the Stirling numbers of the first kind, J. Combin. Theory Ser. A 71 (1995), 343–351.

[14] Z. Kabluchko, V. Vysotsky, and D. Zaporozhets, Convex hulls of random walks, hyperplane arrangements, and Weyl chambers, preprint, 2015,

http://arxiv.org/abs/1510.04073.

[15] Z. Kabluchko, A. Marynych, and H. Sulzbach, Mathematica notebook, 2016,

    http://www.math.uni-muenster.de/statistik/kabluchko/files/stirling.pdf
    and http://www.math.uni-muenster.de/statistik/kabluchko/files//stirling.nb.

[16] Z. Kabluchko, A. Marynych, and H. Sulzbach, General Edgeworth expansions with applications to profiles of random trees, preprint,

    2016,  http://arxiv.org/abs/1606.03920, version with full proofs at
    http://www.math.uni-muenster.de/statistik/kabluchko/files/edgeworth_full.pdf.

[17] L. Kuipers and H. Niederreiter, Uniform Distribution of Sequences, Pure and Applied Mathematics, Wiley-Interscience, 1974.

[18] G. Louchard, Asymptotics of the Stirling numbers of the first kind revisited: a saddle point approach, Discrete Math. Theor. Comput. Sci. 12 (2010), 167–184.

[19] I. Mező, On the maximum of r-Stirling numbers, Adv. in Appl. Math. 41 (2008), 293–306.

[20] L. Moser and M. Wyman, Asymptotic development of the Stirling numbers of the first kind, J. London Math. Soc. 33 (1958), 133–146.

[21] A. Rényi, On the outliers of a series of observations, Magyar Tud. Akad. Mat. Fiz. Oszt. Közl. 12 (1962), 105–121.

[22] M. Sibuya, Log-concavity of Stirling numbers and unimodality of Stirling distributions, Ann. Inst. Statist. Math. 40 (1988), 693–714.

[23] N. M. Temme, Asymptotic estimates of Stirling numbers, Stud. Appl. Math. 89 (1993), 233–243.

[24] H. S. Wilf, The asymptotic behavior of the Stirling numbers of the first kind, J. Combin. Theory Ser. A 64 (1993), 344–349.

[25] H. Yamato, Edgeworth expansions for the number of distinct components associated with the Ewens sampling formula, J. Japan Statist. Soc. 43 (2013), 17–28.

2010 Mathematics Subject Classification: Primary 11B73; Secondary 60C05, 41A60, 60F05, 60F10. Keywords: asymptotic expansion, Ewens distribution, Stirling number of the first kind, mode of a distribution.

(Concerned with sequences \underline{A008275}, \underline{A048994}, \underline{A132393}.)

Received September 30 2016. Revised version received November 21 2016. Published in Journal of Integer Sequences, November 25 2016.

Return to Journal of Integer Sequences home page.