# A refinement of the Ewens sampling formula

Eugene Strahov  Department of Mathematics, The Hebrew University of Jerusalem,
Givat Ram, Jerusalem 91904, Israel<EMAIL_ADDRESS>
###### Abstract.

We consider an infinitely-many neutral allelic model of population genetics
where all alleles are divided into a finite number of classes, and each class
is characterized by its own mutation rate. For this model the allelic
composition of a sample taken from a very large population of genes is
characterized by a random matrix, and the problem is to describe the joint
distribution of the matrix entries. The answer is given by a new
generalization of the classical Ewens sampling formula called the refined
Ewens sampling formula in the present paper. We discuss a Poisson
approximation for the refined Ewens sampling formula, and present its
derivation by several methods. As an application we obtain limit theorems for
the numbers of alleles in different asymptotic regimes.

###### Key words and phrases:

The Ewens sampling formula, the Poisson-Dirichlet distribution and its
generalizations, urn models, genetic population models, random matrices in
population genetics.

##  1\. Introduction

In the context of population genetics the classical Ewens sampling formula [7]
arises as follows. It is assumed that a very large population of genes evolves
according to the neutral Wright-Fisher model with mutations. The population is
haploid (each individual has exactly one parent), and everyone has an equal
chance of reproductive success. The mutations are characterized by a constant
mutation rate, and each mutation leads to a new allele111The paper uses a few
elementary notions from population genetics. In particular, by allele one
means an alternative form of a gene that arises by mutation. For a brief
description of genetic background we refer the reader to Feng [12], Chapter 1,
to Durrett [5], Section 1.1, and to references therein. not previously seen in
the population (this is called the infinite allelic model). If a sample of
size $n$ is taken from such a population, the probability
$p\left(a_{1},\ldots,a_{n};\theta\right)$ that the sample contains $a_{1}$
alleles represented once, $a_{2}$ alleles represented twice, and so on is
given by the Ewens sampling formula

(1.1)
$p\left(a_{1},\ldots,a_{n};\theta\right)=\frac{n!}{(\theta)_{n}}\prod\limits_{j=1}^{n}\left(\frac{\theta}{j}\right)^{a_{j}}\frac{1}{a_{j}!}.$

Here $\theta>0$ is a parameter determined by the mutation rate,
$(\theta)_{n}=\theta(\theta+1)\ldots(\theta+n-1)$ is the Pochhammer symbol,
and $(a_{1},\ldots,a_{n})$ is an allelic partition satisfying
$\sum_{j=1}^{n}ja_{j}=n$. For details of such a derivation of the Ewens
sampling formula we refer the reader to the book by Feng [12], Chapters 2 and
4, to the book by Durrett [5], Section 1.3, and to the book by Etheridge [6],
Chapter 2. For a relation of the Ewens sampling formula to different topics of
mathematics, and for different recent developments around this formula see the
articles by Tavar$\acute{\mbox{e}}$ [20], by Crane [4], and by Griffiths and
Lessard [8].

In the derivation of the Ewens sampling formula the mutations to all possible
alleles are characterized by the same mutation rate $\mu$. In other words, all
possible alleles belong to the same class defined by $\mu$. However, it is not
hard to imagine a situation where all possible alleles can be classified into
finite number of distinct classes. One then assumes that the number of such
classes is $k$, and mutations to alleles from the $l$th class,
$l\in\\{1,\ldots,k\\}$, take place with the mutation rate $\mu_{l}$ associated
with this specific class.

An example motivating to consider infinitely-many allelic models where all
possible alleles are separated into finite number of distinct classes arises
in the study of Cystic Fibrosis (CF). This is a disease which affects
approximately 28000 patients in the United States and approximately 36000
patients in Europe, see O’Reilly and Elphick [17], and references therein. It
is known that the cystic fibrosis is caused by mutations of a specific gene.
This gene identified in 1989 encodes the CF transmembrane conductance
regulator (CFTR) protein, and has more than 1900 mutations with the potential
to cause disease. The mutations can be classified to six different classes,
according to their consequences on CFTR properties, see Table 1. Each class
contains hundreds of different alleles, and characterized by its own mutation
rate.

Class | Consequence of mutation
---|---
1 | No CFTR synthesis
2 | Failure of CFTR to traffic to the cell surface
3 | Failure of CFTR to open at cell surface
4 | Decreased time CFTR channel is open at cell surface
5 | Decreased amount of CFTR at cell surface
6 | Unstable CFTR with decreased half-line at the cell surface

Table 1.  Six different classes of CFTR mutations. The table is taken from
O’Reilly and Elphick [17].

Thus we consider a population of size $2N$ which evolves under the influence
of mutation and random sampling. We assume that each mutation gives rise to a
new allele previously not seen in the population, and all possible alleles
form a countably-infinite set $\mathcal{A}$. The set $\mathcal{A}$ is a union
of $k$ countably-infinite non-intersecting subsets $\mathcal{A}^{(l)}$,
$\mathcal{A}=\cup_{l=1}^{k}\mathcal{A}^{(l)}$. Alleles of class $l$ are the
elements of $\mathcal{A}^{(l)}$, and a mutation to an allele from class $l$ is
characterized by the mutation rate $\mu_{l}$.

Suppose we examine a sample of $n$ gametes taken from such a population. Then
as data we may consider the matrix

(1.2) $\left(\begin{array}[]{ccc}a_{1}^{(1)}&\ldots&a_{1}^{(k)}\\\ \vdots&&\\\
a_{n}^{(1)}&\ldots&a_{n}^{(k)}\end{array}\right),\;\;\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}ja_{j}^{(l)}=n,$

where $a_{j}^{(l)}$ is the number of alleles of class $l$ represented $j$
times in the sample. This data will appear if the experimental techniques
permit a labelling of the alleles we can distinguish according to their
classes only. We say that matrix (1.2) describes the allelic composition of a
sample.

If a random sample of size $n$ is taken from the population, the corresponding
data is a random matrix $A(n)$ of the same format as (1.2). The present paper
initiates a study of such random matrices, and the first problem is to
determine the probability distribution of the matrix entries $A_{j}^{(l)}(n)$
of $A(n)$. As $N\rightarrow\infty$ with $\theta_{1}=4N\mu_{1}$, $\ldots$,
$\theta_{k}=4N\mu_{k}$ fixed, we argue that the joint probability distribution
of $A_{j}^{(l)}(n)$ becomes

(1.3)
$\begin{split}&\operatorname{Prob}\left(A_{j}^{(l)}(n)=a_{j}^{(l)};\;l=1,\ldots,k;\;j=1,\ldots,n\right)\\\
&=\mathbf{1}\left(\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}ja_{j}^{(l)}=n\right)\frac{n!}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\prod\limits_{l=1}^{k}\prod\limits_{j=1}^{n}\left(\frac{\theta_{l}}{j}\right)^{a_{j}^{(l)}}\frac{1}{a_{j}^{(l)}!},\end{split}$

where $\\{a_{j}^{(l)}:\;l=1,\ldots,k;\;j=1,\ldots,n\\}$ are non-negative
integers, and $(\theta)_{n}=\theta(\theta+1)\ldots(\theta+n-1)$ is the
Pochhammer symbol. In equation (1.3) $\theta_{1}>0$, $\ldots$, $\theta_{k}>0$
are parameters determined by the mutation rates $\mu_{1}$, $\ldots$,
$\mu_{k}$. These parameters will be called the mutation parameters throughout
the paper.

For $k=1$ (which means that there is no distinction between alleles from
different classes) formula (1.3) turns into the Ewens sampling formula, and we
refer to (1.1) as to the refined Ewens sampling formula.

The paper is organized as follows. In Section 2 we prove that equation (1.3)
indeed defines a probability distribution on the set of non-negative integers,
see Proposition 2.1. To prove Proposition 2.1 we interpret each matrix (1.2)
as a multiple partition $\Lambda_{n}^{(k)}$ of $n$ into $k$ components, and
equation (1.3) as a measure
$M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$ on the set
$\operatorname{\mathbb{Y}}_{n}^{(k)}$ of all such multiple partitions
$\Lambda_{n}^{(k)}$. Also, we show that expressions (1.3) are consistent for
different values of $n$: if from a sample of size $n$ whose allelic
composition has distribution (1.3) a random sub-sample of size $n-1$ is taken
without replacement, then the distribution of the allelic composition of the
sub-sample is given by (1.3) with $n$ replaced by $n-1$. This property of
(1.3) can be equivalently reformulated as a condition on the system
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$.
Namely,
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
should be a multiple partition structure, see Definition 2.2. We prove this
property of
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
in Proposition 2.3, and thus the consistency of (1.3) for different values of
$n$ follows.

In Section 3 we construct the refined Ewens sampling formula in two different
ways. First, we consider an infinite population divided into $k$ classes,
where each class is a union of a countable number of sub-populations
characterized by their own allelic types. Let $x_{j}^{(l)}$ denote the
proportion of population characterized by the $j$th allelic type from class
$l$. Assuming that $x_{j}^{(l)}$ are fixed numbers we compute the probability
$\mathbb{K}_{n}\left(\Lambda_{n}^{(k)};x^{(1)},\ldots,x^{(k)}\right)$ that the
alleles presented in the sample are characterized by a multiple partition
$\Lambda_{n}^{(k)}$. Then, we allow $x_{j}^{(l)}$ to be random. Under natural
assumptions we obtain that the joint distribution of $x_{j}^{(l)}$ is the
multiple Poisson-Dirichlet distribution
$PD\left(\theta_{1},\ldots,\theta_{k}\right)$. Then the probability that the
alleles presented in the sample are characterized by $\Lambda_{n}^{(k)}$ is
obtained by integration of
$\mathbb{K}_{n}\left(\Lambda_{n}^{(k)};x^{(1)},\ldots,x^{(k)}\right)$ over
certain set $\overline{\nabla}_{0}^{(k)}$, and
$PD\left(\theta_{1},\ldots,\theta_{k}\right)$ is the integration measure on
this set. The result of this integration is the value of
$M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$ at
$\Lambda_{n}^{(k)}$, which is equivalent to the refined Ewens sampling
formula, see Proposition 3.1.

Second, we derive the refined Ewens sampling formula as that for the
stationary distribution (in the infinite population limit) of a Wright-Fisher
type model with mutations. Each mutation leads to a new allele not previously
seen in the population (as in the usual infinite alleles model). The new
assumption here is that all possible alleles are divided into $k$ different
classes, and each class is characterized by its own mutation rate. The model
under considerations leads to the genealogical process which is the coalescent
with killing, the killing is characterized by $k$ different rates (depending
on the allelic class). In the infinite population limit the genealogical
process is associated with a generalization of the Hoppe urn model [9, 10].
The relation with the generalized Hoope urn model enables to compute the
stationary distribution, and thus to derive the refined Ewens sampling
formula, equation (1.3).

In Section 4 we relate the refined Ewens sampling formula with a family of
probability distributions on wreath products. Namely, in Section 4.1 we show
that the refined Ewens sampling formula can be understood as a probability
distribution $P_{t_{1},\ldots,t_{k};n}^{\operatorname{Ewens}}$ on a wreath
product $G\sim S(n)$ of a finite group $G$ with the symmetric group $S(n)$,
see Proposition 4.1. For the wreath products $G\sim S(n)$ we construct an
analogue of the Chinese restaurant process [1] which generates
$P_{t_{1},\ldots,t_{k};n}^{\operatorname{Ewens}}$, see Section 4.2.

In Sections 5 and 6 we present several applications of the refined Ewens
sampling formula. In Section 5 we study the numbers $K_{n}^{(l)}$ of alleles
of class $l$ in a random sample of size $n$. These numbers can be
experimentally observed, and therefore $K_{n}^{(l)}$ are important
probabilistic quantities of interest. We obtain exact formulas for the
expectation of $K_{n}^{(l)}$, for the variance of $K_{n}^{(l)}$, and for the
joint distribution of $K_{n}^{(1)}$, $\ldots$, $K_{n}^{(k)}$. Then we
investigate the asymptotic properties of these random variables as
$n\rightarrow\infty$, in the situation where the mutation parameters are
increasing together with the sample size $n$. Our motivation to consider the
regimes where both $n$ and the mutation parameters tend to infinity comes from
the fact that the mutation parameters are proportional to the population size
assumed very large in the original infinitely-many allelic model. We discuss
convergence in probability and show asymptotic normality of the random
variables $K_{n}^{(l)}$ in different asymptotic regimes, see Theorem 5.2,
Theorem 5.5 and Theorem 5.7.

As it is mentioned above, a random sample can be described by a random matrix
$A(n)$ whose entries are random variables with joint distribution (1.3). In
Section 6 we show that as $n\rightarrow\infty$ the matrix $A(n)$ can be
approximated by a random matrix whose entries are independent Poisson random
variables, see Theorem 6.2. Theorem 6.2 is an extension of the well-known
result by Arratia, Barbour and Tavar$\acute{\mbox{e}}$ [2].

Acknowledgements. I thank Natalia Strahov for bringing to my attention the
paper by O’Reilly and Elphick [17].

##  2\. The refined Ewens sampling formula and multiple partition structures

Formula (1.3) can be understood as a distribution on the set
$\operatorname{\mathbb{Y}}_{n}^{(k)}$ of all multiple partitions
$\Lambda_{n}^{(k)}$ of $n$ into $k$ components. By a multiple partition of $n$
into $k$ components we mean a family
$\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$ of Young
diagrams $\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$ such that the condition

(2.1) $|\lambda^{(1)}|+\ldots+|\lambda^{(k)}|=n$

is satisfied. Here $|\lambda^{(l)}|$ denotes the number of boxes in the Young
diagram $\lambda^{(l)}$. Each matrix (1.2) describing the allelic composition
of a sample of size $n$ can be identified uniquely with a multiple partition
$\Lambda_{n}^{(k)}$ of $n$ into $k$ components. The correspondence is defined
as follows

(2.2) $\left(\begin{array}[]{ccc}a_{1}^{(1)}&\ldots&a_{1}^{(k)}\\\ \vdots&&\\\
a_{n}^{(1)}&\ldots&a_{n}^{(k)}\end{array}\right)\longrightarrow\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right),$

where

(2.3) $\lambda^{(l)}=\left(1^{a_{1}^{(l)}}2^{a_{2}^{(l)}}\ldots
n^{a_{n}^{(l)}}\right),\;l=1,\ldots,k.$

Equation (2.3) says that exactly $a_{j}^{(l)}$ of the rows of $\lambda^{(l)}$
are equal to $j$. It is easy to see that (2.1) is satisfied provided that

(2.4) $\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}ja_{j}^{(l)}=n.$

Thus $\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$ (with
$\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$ defined by equation (2.3)) is an
element of $\operatorname{\mathbb{Y}}_{n}^{(k)}$. Conversely, each
$\Lambda_{n}^{(k)}$ gives rise to a matrix (1.2) describing the allelic
composition of a sample, and we say that $\Lambda_{n}^{(k)}$ defines the
allelic composition of a sample.

######  Proposition 2.1.

Assume that each $\Lambda_{n}^{(k)}\in\operatorname{\mathbb{Y}}_{n}^{(k)}$ is
defined by equations (2.2) and (2.3) in terms of non-negative integers
$a_{j}^{(l)}$ such that the condition (2.4) is satisfied. Set

(2.5)
$M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}(\Lambda_{n}^{(k)})=\frac{n!}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\prod\limits_{l=1}^{k}\prod\limits_{j=1}^{n}\left(\frac{\theta_{l}}{j}\right)^{a_{j}^{(l)}}\frac{1}{a_{j}^{(l)}!}$

We have

$\sum\limits_{\Lambda_{n}^{(k)}\in\operatorname{\mathbb{Y}}_{n}^{(k)}}M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}(\Lambda_{n}^{(k)})=1.$

Thus equation (1.3) indeed defines a probability distribution on the set

$\left\\{a_{j}^{(l)}:\;l=1,\ldots,k;\;j=1,\ldots,n\right\\}$

of non-negative integers.

###### Proof.

From (2.5) we obtain that

(2.6)
$\begin{split}&M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\left(\Lambda_{n}^{(k)}\right)=\frac{\left(\theta_{1}\right)_{|\lambda^{(1)}|}\ldots\left(\theta_{k}\right)_{|\lambda^{(k)}|}}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\\\
&\frac{n!}{|\lambda^{(1)}|!\ldots|\lambda^{(k)}|!}\;M_{\theta_{1};|\lambda^{(1)}|}^{\operatorname{Ewens}}\left(\lambda^{(1)}\right)\ldots
M_{\theta_{k};|\lambda^{(k)}|}^{\operatorname{Ewens}}\left(\lambda^{(k)}\right),\;\;\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right).\end{split}$

Here

(2.7)
$M_{\theta;n}^{\operatorname{Ewens}}\left(\lambda\right)=\frac{n!}{\prod_{j=1}^{n}j^{m_{j}(\lambda)}m_{j}(\lambda)!}\frac{\theta^{l(\lambda)}}{(\theta)_{n}},$

with $\lambda=\left(1^{m_{1}(\lambda)}2^{m_{2}(\lambda)}\ldots\right)$ is the
Ewens probability measure on $\operatorname{\mathbb{Y}}_{n}$, the set of all
Young diagrams with $n$ boxes. Next we use the formula

(2.8)
$n!\sum\limits_{m_{1}+\ldots+m_{k}=n}\frac{(\theta_{1})_{m_{1}}\ldots(\theta_{k})_{m_{k}}}{m_{1}!\ldots
m_{k}!}=\left(\theta_{1}+\ldots+\theta_{k}\right)_{n},$

and the fact that each $M_{\theta^{(l)},m_{l}}^{\operatorname{Ewens}}$ is a
probability measure on $\operatorname{\mathbb{Y}}_{m_{l}}$. ∎

A random multiple partition of $n$ into $k$ components is a random variable
$\Lambda_{n}^{(k)}$ taking values in the set
$\operatorname{\mathbb{Y}}_{n}^{(k)}$. If a sample of $n$ gametes is chosen at
random from a population, then the allelic composition of a sample is
described by the random variable $\Lambda_{n}^{(k)}$.

######  Definition 2.2.

A multiple partition structure222Multiple partition structures were first
introduced in Strahov [18]. Multiple partition structures can be viewed as
generalizations of Kingman partition structures [13, 14].  is a sequence
$\mathcal{M}_{1}^{(k)}$, $\mathcal{M}^{(n)}_{2}$, $\ldots$ of distributions
for $\Lambda_{1}^{(k)}$, $\Lambda_{2}^{(k)}$, $\ldots$ which is consistent in
the following sense: if the allelic composition of a sample of $n$ gametes is
$\Lambda_{n}^{(k)}$, and a random sub-sample of size $n-1$ is taken
independently of $\Lambda_{n}^{(k)}$, then the allelic composition
$\Lambda_{n-1}^{(k)}$ is distributed according to $\mathcal{M}_{n-1}^{(k)}$.

######  Proposition 2.3.

The system
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
is a multiple partition structure.

###### Proof.

It was proved in Strahov [18], Section 8.1 that a sequence
$\mathcal{M}_{1}^{(k)}$, $\mathcal{M}^{(n)}_{2}$, $\ldots$ of distributions
for $\Lambda_{1}^{(k)}$, $\Lambda_{2}^{(k)}$, $\ldots$ is a multiple partition
structure if and only if the consistency condition

(2.9)
$\mathcal{M}_{n-1}^{(k)}\left(\Lambda_{n-1}^{(k)}\right)=\sum\limits_{\Lambda_{n}^{(k)}\in\operatorname{\mathbb{Y}}_{n}^{(k)}}\operatorname{Prob}\left(\Lambda_{n-1}^{(k)}\biggl{|}\Lambda_{n}^{(k)}\right)\mathcal{M}_{n}^{(k)}\left(\Lambda_{n}^{(k)}\right),$

is satisfied for each $\Lambda_{n-1}^{(k)}$,
$\Lambda_{n-1}^{(k)}\in\operatorname{\mathbb{Y}}_{n-1}^{(k)}$. Here the
transition probability
$\operatorname{Prob}\left(\Lambda_{n-1}^{(k)}\biggl{|}\Lambda_{n}^{(k)}\right)$
to get $\Lambda_{n-1}^{(k)}=\left(\mu^{(1)},\ldots,\mu^{(k)}\right)$ from
$\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$ is defined
by

(2.10)
$\operatorname{Prob}\left(\Lambda_{n-1}^{(k)}\biggl{|}\Lambda_{n}^{(k)}\right)=\left\\{\begin{array}[]{ll}\frac{1}{n}m_{L^{(1)}}(\lambda^{(1)})L^{(1)},&\lambda^{(1)}\searrow\mu^{(1)},\\\
\frac{1}{n}m_{L^{(2)}}(\lambda^{(2)})L^{(2)},&\lambda^{(2)}\searrow\mu^{(2)},\\\
\vdots,&\\\
\frac{1}{n}m_{L^{(k)}}(\lambda^{(k)})L^{(k)},&\lambda^{(k)}\searrow\mu^{(k)},\\\
0,&\mbox{otherwise},\end{array}\right.$

where the notation $\lambda^{(l)}\searrow\mu^{(l)}$ means that $\mu^{(l)}$ is
obtained from $\lambda^{(l)}$ by extracting one box, $L^{(l)}$ is the length
of the row in $\lambda^{(l)}$ from which a box is removed to get $\mu^{(l)}$,
and $m_{L^{(l)}}\left(\lambda^{(l)}\right)$ is the number of rows of size
$L^{(l)}$ in $\lambda^{(l)}$.

In order to see that the system
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
is a multiple partition structure we first check directly using (2.7) that

(2.11)
$M_{\theta,n-1}^{\operatorname{Ewens}}(\mu)=\sum\limits_{\lambda\searrow\mu}\left(\frac{1}{n}m_{L}(\lambda)L\right)M_{\theta,n}^{\operatorname{Ewens}}(\lambda).$

The fact that the system
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
satisfies (2.9) follows from equations (2.6), (2.10) and (2.11) by direct
calculations. ∎

Thus
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
is a multiple partition structure in the sense of Definition 2.2. In what
follows we refer to
$\left(M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\right)_{n=1}^{\infty}$
as to the Ewens multiple partition structure.

##  3\. A construction of the Ewens multiple partition structure

###  3.1. The Ewens multiple partition structure and the multiple Poisson-
Dirichlet distribution

Consider an infinite population divided into $k$ classes, where each class is
a union of a countable number of sub-populations characterized by their own
allelic types. Thus a sub-population is labeled by the pair $(j,l)$ of
numbers, indicating that it is characterized by $j$th allelic type from the
mutation class numbered by $l$. Here $j=1,2,\ldots,$ and $l=1,\ldots,k$.

Suppose that the proportion of the population labeled by $(j,l)$ is
$x_{j}^{(l)}$. Then $x_{j}^{(l)}\geq 0$, and

(3.1) $\begin{split}&x_{1}^{(1)}+x_{2}^{(1)}+\ldots=\delta^{(1)}\\\ \vdots\\\
&x_{1}^{(k)}+x_{2}^{(k)}+\ldots=\delta^{(k)}.\end{split}$

The parameter $\delta^{(l)}$ characterizes the frequency of the $l$th class in
the population, so the condition

(3.2) $\delta^{(1)}+\ldots+\delta^{(k)}=1$

should be satisfied.

Let $n$ be the size of a sample chosen from the population. The alleles
presented in the sample can be characterized by a multiple partition

(3.3)
$\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right),\;\;|\lambda^{(1)}|+\ldots+|\lambda^{(k)}|=n.$

Recall that the meaning of $\Lambda_{n}^{(k)}$ is the following. Write

(3.4) $\lambda^{(l)}=\left(1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots
n^{m_{n}^{(l)}}\right),$

for each $l\in\left\\{1,\ldots,k\right\\}$. Then $m_{j}^{(l)}$ is the number
of alleles of the $l$th class that appear $j$ times in the sample. Here we
explain that the distribution of $\Lambda_{n}^{(k)}$ should be the multiple
Ewens distribution.

Assume first that the frequencies $x_{j}^{(l)}$ introduced above are fixed
positive numbers satisfying (3.1). The probability that for each
$l=1,\ldots,k$ a random sample of size $n$ has $m_{1}^{(l)}$ of alleles of the
class $l$ represented once, $m_{2}^{(l)}$ alleles of the class $l$ represented
twice and so on is

(3.5)
$\begin{split}&\mathbb{K}\left(\Lambda_{n}^{(k)},x^{(1)},\ldots,x^{(k)}\right)=\frac{n!}{\prod_{l=1}^{k}\prod_{j=1}^{n}\left(j!\right)^{m_{j}^{(l)}}}\\\
&\times\sum\left[\left(x_{1}^{(1)}\right)^{\nu_{1}^{(1)}}\left(x_{2}^{(1)}\right)^{\nu_{2}^{(1)}}\ldots\right]\ldots\left[\left(x_{1}^{(k)}\right)^{\nu_{1}^{(k)}}\left(x_{2}^{(k)}\right)^{\nu_{2}^{(k)}}\ldots\right],\end{split}$

where the sum is over all $k$ sequences
$\left(\nu_{1}^{(1)},\nu_{2}^{(1)},\ldots\right)$, $\ldots$,
$\left(\nu_{1}^{(k)},\nu_{2}^{(k)},\ldots\right)$ of integers
$0\leq\nu_{j}^{(l)}\leq n$ such that the number of elements equal to $j$ in
$\left(\nu_{1}^{(l)},\nu_{2}^{(l)},\ldots\right)$ is precisely $m_{j}^{(l)}$,
where $l=1,\ldots,k$, and $j=1,\ldots,n$. We can write

(3.6)
$\begin{split}&\mathbb{K}\left(\Lambda_{n}^{(k)},x^{(1)},\ldots,x^{(k)}\right)=\frac{n!}{\prod_{l=1}^{k}\prod_{j=1}^{n}\left(j!\right)^{m_{j}^{(l)}}}\;m_{\lambda^{(1)}}\left(x^{(1)}\right)\ldots
m_{\lambda^{(k)}}\left(x^{(k)}\right),\end{split}$

where $m_{\lambda^{(l)}}\left(x^{(l)}\right)$ is the monomial symmetric
function, and $\Lambda_{n}^{(k)}$ is defined by (3.3) and (3.4).

Next, we allow $x^{(1)}$, $\ldots$, $x^{(k)}$ to be random. If we restrict
ourselves to samples containing alleles from class $l$, then the corresponding
frequencies should by governed by the Poisson-Dirichlet distribution333For the
definition of the Poisson-Dirichlet distribution, and its connection to the
Ewens sampling formula we refer the reader to book Feng [12].
$PD(\theta_{l})$, where $\theta_{l}$ is the mutation rate for the class $l$.
The parameters $\delta_{1}$, $\ldots$, $\delta_{k}$ are random variables, and
are interpreted as proportions of sub-populations. Clearly, these random
variables should satisfy

(3.7)
$\mathbb{E}\left(\delta_{1}\right)=\frac{\theta_{1}}{\theta_{1}+\ldots+\theta_{k}},\ldots,\mathbb{E}\left(\delta_{k}\right)=\frac{\theta_{k}}{\theta_{1}+\ldots+\theta_{k}}.$

The Dirichlet distribution $D\left(\theta_{1},\ldots,\theta_{k}\right)$ is a
suitable choice for the joint distribution of the random variables
$\delta_{1}$, $\ldots$, $\delta_{k}$.

Therefore, it is reasonable to suggest the following description of the random
frequencies $x^{(l)}=\left(x_{1}^{(l)},x_{2}^{(l)},\ldots\right)$,
$l=1,\ldots,k$. Let us assume that for each $l=1,\ldots,k$ the sequence
$y^{(l)}=\left(y_{1}^{(l)}\geq y_{2}^{(l)}\ldots\geq 0\right)$ of random
variables is governed by the Poisson-Dirichlet distribution by
$PD\left(\theta_{l}\right)$, and the sequences $y^{(1)}$, $\ldots$, $y^{(k)}$
are independent. In addition, let $\delta_{1}$, $\ldots$, $\delta_{k}$ be
random variables independent on $y^{(1)}$, $\ldots$, $y^{(k)}$ whose joint
distribution is $D(\theta_{1},\ldots,\theta_{k})$. Then each sequence
$x^{(l)}$ is distributed as

$\delta_{l}y^{(l)}=\left(\delta_{l}y_{1}^{(l)}\geq\delta_{l}y_{2}^{(l)}\geq\ldots\geq
0\right).$

The joint distribution of $x^{(1)}$, $\ldots$, $x^{(k)}$ is called the
multiple Poisson-Dirichlet distribution, and it is denoted by
$PD\left(\theta_{1},\ldots,\theta_{k}\right)$.

######  Proposition 3.1.

We have

(3.8)
$\int\limits_{\overline{\nabla}_{0}^{(k)}}\mathbb{K}\left(\Lambda_{n}^{(k)},\omega\right)PD(\theta_{1},\ldots,\theta_{k})(d\omega)=M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\left(\Lambda_{n}^{(k)}\right),$

where $\overline{\nabla}_{0}^{(k)}$ is defined by

(3.9)
$\begin{split}\overline{\nabla}^{(k)}_{0}=&\biggl{\\{}(x,\delta)\biggr{|}x=\left(x^{(1)},\ldots,x^{(k)}\right),\;\delta=\left(\delta^{(1)},\ldots,\delta^{(k)}\right);\\\
&x^{(l)}=\left(x^{(l)}_{1}\geq x^{(l)}_{2}\geq\ldots\geq
0\right),\;\delta^{(l)}\geq 0,\;1\leq l\leq k,\\\
&\mbox{where}\;\;\sum\limits_{i=1}^{\infty}x_{i}^{(l)}=\delta^{(l)},\;1\leq
l\leq
k,\;\mbox{and}\;\sum\limits_{l=1}^{k}\delta^{(l)}=1\biggl{\\}}.\end{split}$

###### Proof.

The integration variables in the left-hand side of equation (3.8) can be
represented as $x^{(l)}=\delta_{l}y^{(l)}$, $l=1,\ldots,k$. The joint
distribution of $\delta_{1}$, $\ldots$, $\delta_{k}$ is
$D\left(\theta_{1},\ldots,\theta_{k}\right)$, and the distribution of
$y^{(l)}$ is $PD(\theta_{l})$. The left-hand side of equation (R1.7) can be
written as

(3.10)
$\begin{split}\frac{n!}{\prod_{l=1}^{k}\prod_{j}^{n}(j!)^{m_{j}^{(l)}}}\int\limits_{\overline{\nabla}_{0}^{(1)}}\ldots\int\limits_{\overline{\nabla}_{0}^{(1)}}&\underset{\delta_{1}+\ldots+\delta_{k}=1}{\underset{\delta_{1}\geq
0,\ldots,\delta_{k}\geq
0}{\int\ldots\int}}\left(\prod\limits_{l=1}^{k}m_{\lambda^{(l)}}\left(\delta_{l}y^{(l)}\right)PD\left(\theta_{l}\right)(dy^{(l)})\right)\\\
&\times D(\theta_{1},\ldots,\theta_{k})\left(d\delta_{1}\ldots
d\delta_{k}\right).\end{split}$

Note that

(3.11)
$m_{\lambda^{(l)}}\left(\delta_{l}y^{(l)}\right)=\left(\delta_{l}\right)^{|\lambda^{(l)}|}m_{\lambda^{(l)}}\left(y^{(l)}\right),$

and that

(3.12)
$\begin{split}&\underset{\delta_{1}+\ldots+\delta_{k}=1}{\underset{\delta_{1}\geq
0,\ldots,\delta_{k}\geq
0}{\int\ldots\int}}\delta_{1}^{|\lambda^{(1)}|}\ldots\delta_{k}^{|\lambda^{(k)}|}D(\theta_{1},\ldots,\theta_{k})\left(d\delta_{1}\ldots
d\delta_{k}\right)=\frac{(\theta_{1})_{|\lambda^{(1)}|}\ldots(\theta_{k})_{|\lambda^{(k)}|}}{(\theta_{1}+\ldots+\theta_{k})_{n}}.\end{split}$

Also, the integral representation of the Ewens distribution
$M_{\theta_{l},n}^{\operatorname{Ewens}}$ is

(3.13)
$\frac{|\lambda^{(l)}|!}{\prod_{j=1}^{n}(j!)^{m_{j}^{(l)}}}\int\limits_{\overline{\nabla}_{0}^{(1)}}m_{\lambda^{(l)}}\left(y^{(l)}\right)PD(\theta_{l})(dy^{(l)})=M_{\theta_{l},|\lambda^{(i)}|}^{\operatorname{Ewens}}\left(\lambda^{(l)}\right),$

see Kingman [13]. Therefore, the left-hand side of equation (3.8) is equal to

$\frac{(\theta_{1})_{|\lambda^{(1)}|}\ldots(\theta_{k})_{|\lambda^{(k)}|}}{(\theta_{1}+\ldots+\theta_{k})_{n}}\frac{n!}{|\lambda^{(1)}|!\ldots|\lambda^{(k)}|!}M_{\theta_{1},|\lambda^{(1)}|}^{\operatorname{Ewens}}\left(\lambda^{(1)}\right)\ldots
M_{\theta_{k},|\lambda^{(k)}|}^{\operatorname{Ewens}}\left(\lambda^{(k)}\right),$

which is
$M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\left(\Lambda_{n}^{(k)}\right)$,
$\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$. ∎

Note that equation (3.8) provides an integral representation for the Ewens
multiple partition structure
$\left(M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n}\right)_{n=1}^{\infty}$.

###  3.2. The coalescent with killing and a generalization of the Hoppe urn
model

In this Section the refined Ewens sampling formula is obtained as a stationary
distribution for a model where a population evolves under the influence of
mutation and random sampling. Our arguments here are generalizations of those
from Feng [12], Chapter 4, and from Durrett [5], Section 1.3.

To begin with, let us assume that the size of the population is $2N$, that the
population is habloid (i.e. each individual has just one parent), and that the
population evolves according to the Wright-Fisher model with no mutation. If a
sample of size $n$ is taken from this population, denote by $P_{nm}$ the
probability that the number of the ancestors of this sample is $m$. There is a
simple formula for $P_{nm}$, namely

(3.14) $P_{nm}=S(n,m)\frac{(2N)(2N-1)\ldots(2N-m+1)}{(2N)^{n}},$

where $S(n,m)$ is the total number of ways of partitioning $\\{1,\ldots,n\\}$
into $m$ non-empty subsets. We have $S(n,n)=1$, and
$S(n,n-1)=\left(\begin{array}[]{c}n\\\ 2\end{array}\right)$. Note that

$P_{nn}=1-\frac{\left(\begin{array}[]{c}n\\\
2\end{array}\right)}{2N}+O\left(\frac{1}{(2N)^{2}}\right),\;\;P_{nn-1}=\frac{\left(\begin{array}[]{c}n\\\
2\end{array}\right)}{2N}+O\left(\frac{1}{(2N)^{2}}\right),$

and

$P_{nm}=O\left(\frac{1}{(2N)^{2}}\right),\;\;m=1,2,\ldots,n-2.$

Next, we say that mutations are possible, and that each mutation leads to a
new allele not previously seen in the population. All possible alleles are
divided into $k$ different classes characterized by their own mutation rates
$\mu_{1}$, $\ldots$, $\mu_{k}$. Each $\mu_{l}$, $l=1,\ldots,k$, is understood
as the probability per individual to mutate into a new allele of class $l$,
$l\in\left\\{1,\ldots,k\right\\}$. The mutation rates are assumed to be of
order $O\left(\frac{1}{N}\right)$, so we write

$\mu_{1}=\frac{\theta_{1}}{4N},\ldots,\mu_{k}=\frac{\theta_{k}}{4N},$

where $\theta_{1}>0$, $\ldots$, $\theta_{k}>0$ are fixed numbers. These
considerations lead to a Markov chain $Z_{m}$ defined as follows. Set
$Z_{0}=n$ (recall that $n$ is the size of the sample under considerations),
$Z_{1}$ is the total number of ancestors, one generation back, and $Z_{m}$ is
the total number of ancestors, $m$ generations back. The transition
probability $\mathcal{P}_{nm}$ of this Markov chain is defined by

(3.15)
$\mathcal{P}_{nm}=\left\\{\begin{array}[]{ll}1-\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{4N},&m=n,\\\
\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{4N},&m=n-1\\\
0,&\hbox{otherwise.}\end{array}\right.$

In particular, the transition from the state $n$ to the state $n-1$ can occur
due to the following reasons. The first reason is that two individuals from
the generation of size $n$ will pick the same parent, the probability of this
event is $\left(\begin{array}[]{c}n\\\ 2\end{array}\right)\frac{1}{2N}$ (where
$\left(\begin{array}[]{c}n\\\ 2\end{array}\right)$ is the number of ways of
picking two of the $n$ individuals, and $\frac{1}{2N}$ is the probability that
two individuals have the same parent). The second reason is the annihilation:
the number $\frac{\theta_{l}n}{4N}$ gives the probability of the annihilation
of an allele of class $l$, $l\in\left\\{1,\ldots,k\right\\}$.

Denote by $\tau_{n}$ the number of steps the Markov chain spends at the state
$n$. Then

(3.16)
$\operatorname{Prob}\left\\{\tau_{n}=m\right\\}=\left[1-\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{4N}\right]^{m-1}\left[\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{4N}\right].$

This gives

(3.17) $\begin{split}&\operatorname{Prob}\left\\{\tau_{n}\leq
2Nt\right\\}=\left[\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{4N}\right]\sum\limits_{m=1}^{[2Nt]}\left[1-\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{4N}\right]^{m-1}.\end{split}$

It follows that

(3.18)
$\underset{N\rightarrow\infty}{\lim}\operatorname{Prob}\left\\{\tau_{n}\leq
2Nt\right\\}=1-\exp{\left[-\left(\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{2}\right)t\right]}.$

We conclude that $\frac{\tau_{n}}{2N}$ turns into the exponential random
variable $T_{n}$ with rate

$\frac{n(n-1)+\theta_{1}n+\ldots+\theta_{k}n}{2},$

as $N\rightarrow\infty$.

The Markov chain described above can be approximated by a pure-death process
$D_{t}$, where $D_{t}$ denotes the number of ancestors at time $t$. The
process starts at $n$, and we have

$\operatorname{Prob}\left(D_{t+h}=n-1\biggl{|}D_{t}=n\right)=\frac{n(n-1)}{2}+\frac{\theta_{1}n}{2}+\ldots+\frac{\theta_{k}n}{2}+o(h),$

as $h\rightarrow 0$. The expression

$\frac{n(n-1)}{2}+\frac{\theta_{1}n}{2}+\ldots+\frac{\theta_{k}n}{2}$

is the transition intensity of the process $D_{t}$. It is interpreted as we
have a coalescing event with probability

(3.19) $\frac{n-1}{n-1+\theta_{1}+\ldots+\theta_{k}},$

or a mutation into an allele of class $l$ with probability

(3.20) $\frac{\theta_{l}}{n-1+\theta_{1}+\ldots+\theta_{k}},$

where $l=1,\ldots,k$. The coalescent with mutation corresponds to a
generalized Hoppe’s urn model described below.

A generalization of the Hoppe urn model. Consider an urn which contains
objects of $k$ different classes (for example, balls, cubs, pyramids, and so
on). For each $l\in\\{1,\ldots,k\\}$ the $l$th class is represented by a black
object with mass $\theta_{l}>0$, and by a number of colored objects of the
same class. We agree that every colored (non-black) object has mass 1.

Assume that such an urn is given, and it contains $n$ objects. A multiple
partition $\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$
of $n$ into $k$ components can be associated to the urn as follows. Consider
the colored objects of the first class, and denote by $\lambda_{1}^{(1)}$ the
number of such objects of the most frequent color, by $\lambda_{2}^{(1)}$ the
number of such objects of the next frequent color, and so on. As a result, we
will obtain a Young diagram $\lambda^{(1)}$ associated with the objects of the
first class in the urn. The Young diagrams $\lambda^{(2)}$, $\ldots$,
$\lambda^{(k)}$ associated with the objects of other classes are constructed
in the same way.

At each time, an object is selected at random with a probability proportional
to its mass. If a colored (non-black) object is selected, that object and
another object of the same class and of the same color are returned to the
urn. If a black object is chosen, that object is returned to the urn with an
object of the same class, and of a new color. This new object has mass $1$.

Let $\Lambda_{n}^{(k)}$ be a multiple partition associated with the urn
containing objects of $k$ different classes (as described above). Then
$\Lambda_{n}^{(k)}$ is a random variable, and the collection
$\left\\{\Lambda_{n}^{(k)}\right\\}_{n=1}^{\infty}$ is a Markov process. The
marginal distribution of this process is given by the next Theorem.

######  Theorem 3.2.

Let $\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$ be fixed Young diagrams such
that

$\left|\lambda^{(1)}\right|+\ldots+\left|\lambda^{(k)}\right|=n.$

Then $\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$ is a fixed multiple
partition of $n$ into $k$ components. We have

(3.21)
$\operatorname{Prob}\left(\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)\right)=M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right),$

where the probability measure
$M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$ is defined by
equations (2.6) and (2.7).

###### Proof.

We proceed by induction. If $n=1$, the possible multiple partitions are

$((1),\varnothing,\ldots,\varnothing),\;(\varnothing,(1),\varnothing,\ldots,\varnothing),\ldots,(\varnothing,\varnothing,\ldots,(1)).$

In particular, the multiple partition

$\left(\underset{l}{\underbrace{\varnothing,\ldots,\varnothing,(1)}},\varnothing,\ldots,\varnothing\right)$

has probability $\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}}$. On the
other hand,

(3.22)
$\begin{split}&M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n=1}\left(\left(\underset{l}{\underbrace{\varnothing,\ldots,\varnothing,(1)}},\varnothing,\ldots,\varnothing\right)\right)\\\
&=\frac{(\theta_{l})_{1}}{(\theta_{1}+\ldots+\theta_{k})_{1}}\frac{1!}{0!\ldots
1!\ldots
0!}M^{\operatorname{Ewens}}_{\theta_{1},|\lambda^{(1)}|=0}\left(\varnothing\right)\ldots
M^{\operatorname{Ewens}}_{\theta_{l},|\lambda^{(l)}|=1}\left((1)\right)\ldots
M^{\operatorname{Ewens}}_{\theta_{k},|\lambda^{(k)}|=0}\left(\varnothing\right)\\\
&=\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}},\end{split}$

so the result holds true.

Denote by $P_{n}^{\operatorname{Hoppe}}$ the probability measure describing
the distribution of the multiple partition $\Lambda_{n}^{(k)}$. We have

(3.23)
$P_{n}^{\operatorname{Hoppe}}\left(\Lambda_{n}^{(k)}\right)=\sum\limits_{\Lambda_{n-1}^{(k)}\in\operatorname{\mathbb{Y}}_{n-1}^{(k)}}P^{\operatorname{Hoppe}}_{n-1}\left(\Lambda_{n-1}^{(k)}\right)p\left(\Lambda_{n-1}^{(k)},\Lambda_{n}^{(k)}\right),$

where $p\left(\Lambda_{n-1}^{(k)},\Lambda_{n}^{(k)}\right)$ is the transition
probability of the Markov process
$\left\\{\Lambda_{n}^{(k)}\right\\}_{n=1}^{\infty}$. Let us find this
transition probability explicitly. We rewrite equation (3.23) as

(3.24)
$\begin{split}&P^{\operatorname{Hoppe}}_{n}\left(\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)\right)\\\
&=\sum\limits_{\mu^{(1)}\nearrow\lambda^{(1)}}P_{n-1}^{\operatorname{Hoppe}}\left(\left(\mu^{(1)},\lambda^{(2)},\ldots,\lambda^{(k)}\right)\right)p\left(\left(\mu^{(1)},\lambda^{(2)},\ldots,\lambda^{(k)}\right),\left(\lambda^{(1)},\lambda^{(2)},\ldots,\lambda^{(k)}\right)\right)\\\
&+\sum\limits_{\mu^{(2)}\nearrow\lambda^{(2)}}P_{n-1}^{\operatorname{Hoppe}}\left(\left(\lambda^{(1)},\mu^{(2)},\ldots,\lambda^{(k)}\right)\right)p\left(\left(\lambda^{(1)},\mu^{(2)},\ldots,\lambda^{(k)}\right),\left(\lambda^{(1)},\lambda^{(2)},\ldots,\lambda^{(k)}\right)\right)\\\
&+\ldots\\\
&+\sum\limits_{\mu^{(k)}\nearrow\lambda^{(k)}}P_{n-1}^{\operatorname{Hoppe}}\left(\left(\lambda^{(1)},\lambda^{(2)},\ldots,\mu^{(k)}\right)\right)p\left(\left(\lambda^{(1)},\lambda^{(2)},\ldots,\mu^{(k)}\right),\left(\lambda^{(1)},\lambda^{(2)},\ldots,\lambda^{(k)}\right)\right),\end{split}$

where $\mu^{(l)}\nearrow\lambda^{(l)}$ means that $\lambda^{(l)}$ is obtained
from $\mu^{(l)}$ by adding one box. In order to get a formula for the
transition probability

$p\left(\left(\lambda^{(1)},\ldots,\mu^{(l)},\ldots,\lambda^{(k)}\right),\left(\lambda^{(1)},\ldots,\lambda^{(l)},\ldots,\lambda^{(k)}\right)\right)$

(where $\mu^{(l)}\nearrow\lambda^{(l)}$ ) we observe that $\lambda^{(l)}$ can
be obtained from $\mu^{(l)}$ in two ways:
(a) A new box is attached to the bottom of $\mu^{(l)}$. If
$\lambda^{(l)}=\left(1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots\right)$, then
$\mu^{(l)}=\left(1^{m_{1}^{(l)}-1}2^{m_{2}^{(l)}}\ldots\right)$. This
corresponds to the addition of an object of class $l$ with a new color. In
this case, the transition probability for the generalized Hoppe urn model is

$\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+n-1}.$

(b) A row of length $j$ in $\mu^{(l)}$ becomes a row of length $j+1$ in
$\lambda^{(l)}$, where $j$ takes values in the set
$\left\\{1,2,\ldots,\lambda_{1}^{(l)}-1\right\\}$. If

$\lambda^{(l)}=\left(1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}...\right),$

then

$\mu^{(l)}=\left(1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots(j-1)^{m_{j-1}^{(l)}}j^{m_{j}^{(l)}+1}(j+1)^{m_{j+1}^{(l)}-1}(j+2)^{m_{j+2}^{(l)}}\ldots\right).$

In this case, at time $n-1$ a colored (non-black) object of class $l$ is
chosen from the existing $j$ objects of class $l$, and of the same color. The
chosen object is added with an additional object of the same class, and of the
same color. The corresponding transition probability is

$\frac{j\left(m_{j}^{(l)}+1\right)}{\theta_{1}+\ldots+\theta_{k}+n-1}.$

It follows that the distribution $P^{\operatorname{Hoppe}}_{n}$ of the
generalized Hoppe urn satisfies the recursion

$\begin{split}&1=\sum\limits_{l=1}^{k}\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+n-1}\frac{P^{\operatorname{Hoppe}}_{n-1}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;1^{m_{1}^{(l)}-1}2^{m_{2}^{(l)}}\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}{P^{\operatorname{Hoppe}}_{n}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}\\\
&+\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{\lambda_{1}^{(l)}-1}\frac{j(m_{j}^{(l)}+1)P^{\operatorname{Hoppe}}_{n-1}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;\ldots;1^{m_{1}^{(l)}}\ldots
j^{m_{j}^{(l)}+1}(j+1)^{m_{j+1}^{(l)}-1}\ldots;\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}{(\theta_{1}+\ldots+\theta_{k}+n-1)P^{\operatorname{Hoppe}}_{n}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;\ldots;1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots;\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}\end{split}$

Let us consider the probability measure
$M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n}$ defined by
equations (2.6) and (2.7). It is enough to show that
$M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n}$ satisfies the same
recursion as $P^{Hoppe}_{n}$. We have

(3.25)
$\frac{M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n-1}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;1^{m_{1}^{(l)}-1}2^{m_{2}^{(l)}}\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}{M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}=\frac{m_{1}^{(l)}\left(\theta_{1}+\ldots+\theta_{k}+n-1\right)}{n\theta_{l}},$

and

(3.26)
$\begin{split}&\frac{M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n-1}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;\ldots;1^{m_{1}^{(l)}}\ldots
j^{m_{j}^{(l)}+1}(j+1)^{m_{j+1}^{(l)}-1}\ldots;\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}{M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n}\left(1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots;\ldots;1^{m_{1}^{(l)}}2^{m_{2}^{(l)}}\ldots;\ldots;1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots\right)}\\\
&=\frac{\theta_{1}+\ldots+\theta_{k}+n-1}{n}\;\frac{j+1}{j}\frac{m_{j+1}^{(l)}}{m_{j}^{(l)}+1}.\end{split}$

We see that $M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k};n}$
satisfies the same recurrence relation as $P^{\operatorname{Hoppe}}_{n}$
provided that the condition

(3.27)
$1=\sum\limits_{l=1}^{k}\frac{m_{1}^{(l)}}{n}+\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{\lambda_{1}^{(l)}-1}\frac{(j+1)m_{j+1}^{(l)}}{n}.$

is satisfied. Since

$m_{1}^{(l)}+\sum\limits_{j=1}^{\lambda_{1}^{(l)}-1}(j+1)m_{j+1}^{(l)}=\sum\limits_{j=1}^{\lambda_{1}^{(l)}}jm_{j}^{(l)}=|\lambda^{(l)}|,$

and $|\lambda^{(1)}|+\ldots+|\lambda^{(k)}|=n$, condition (3.27) holds true
indeed. ∎

Finally, interpret $\Lambda_{n}^{(k)}$ from Theorem 3.2 as a multiple
partition describing the allelic composition of a sample of size $n$ (see the
beginning of Section 2) whose genealogy is described by coalescence with
probability (3.19), and mutations with probabilities (3.20). With such an
interpretation, Theorem 3.2 says that the distribution of $\Lambda_{n}^{(k)}$
is given by the refined Ewens sampling formula.

###  3.3. Consistency of
$M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$: the reduction to
the Ewens distribution on Young diagrams

The interpretation of $\theta_{1}$, $\ldots$, $\theta_{k}$ as parameters
describing mutations into alleles of different classes dictates a consistency
condition on $M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$.
Assume that $M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$
defines the distribution of $\Lambda_{n}^{(k)}$ describing the allelic
composition of a sample of size $n$. If the experimental techniques does not
enable to distinguish between $k$ different classes of alleles, then the
allelic composition of a sample should be described by a Young diagram with
$n$ boxes distributed according to the usual Ewens formula with the parameter
$\theta_{1}+\ldots+\theta_{k}$. More precisely, the following Proposition
should be valid.

######  Proposition 3.3.

Let $\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$ be Young diagrams. Write
$\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$ as

(3.28) $\begin{split}&\lambda^{(1)}=1^{m_{1}^{(1)}}2^{m_{2}^{(1)}}\ldots
n^{m_{n}^{(1)}},\\\ &\vdots\\\
&\lambda^{(k)}=1^{m_{1}^{(k)}}2^{m_{2}^{(k)}}\ldots
n^{m_{n}^{(k)}},\end{split}$

where $n=|\lambda^{(1)}|+\ldots+|\lambda^{(k)}|$. The union
$\lambda^{(1)}\cup\ldots\cup\lambda^{(k)}$ of $\lambda^{(1)}$, $\ldots$,
$\lambda^{(k)}$ is defined by

$\lambda^{(1)}\cup\ldots\cup\lambda^{(k)}=1^{m_{1}^{(1)}+\ldots+m_{1}^{(k)}}2^{m_{2}^{(1)}+\ldots+m_{2}^{(k)}}\ldots
n^{m_{n}^{(1)}+\ldots+m_{n}^{(k)}}.$

In words, the Young diagram $\lambda^{(1)}\cup\ldots\cup\lambda^{(k)}$ is
obtained from $\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$ by combining together
all rows of $\lambda^{(1)}$, $\ldots$, $\lambda^{(k)}$.

Let $M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}$ be the
probability measure on $\operatorname{\mathbb{Y}}_{n}^{(k)}$ defined by
equation (2.6). Then we have

(3.29)
$\underset{\lambda^{(1)}\cup\ldots\cup\lambda^{(k)}=\lambda}{\sum\limits_{\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)}}M_{\theta_{1},\ldots,\theta_{k};n}^{\operatorname{Ewens}}\left(\Lambda_{n}^{(k)}\right)=M_{\theta_{1}+\ldots+\theta_{k};n}^{\operatorname{Ewens}}\left(\lambda\right),$

where $M_{\theta;n}^{\operatorname{Ewens}}\left(\lambda\right)$ is defined by
equation (2.7).

###### Proof.

Writing the left-hand side of equation (3.29) explicitly we obtain

$\begin{split}&\frac{n!}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\underset{m_{n}^{(1)}+\ldots+m_{n}^{(k)}=m_{n}}{\underset{\vdots}{\sum\limits_{m_{1}^{(1)}+\ldots+m_{1}^{(k)}=m_{1}}}}\frac{\theta_{1}^{m_{1}^{(1)}+\ldots+m_{n}^{(1)}}\ldots\theta_{k}^{m_{1}^{(k)}+\ldots+m_{n}^{(k)}}}{\left(\prod\limits_{j=1}^{n}j^{m_{j}^{(1)}}m_{j}^{(1)}!\right)\ldots\left(\prod\limits_{j=1}^{n}j^{m_{j}^{(k)}}m_{j}^{(k)}!\right)}\\\
&=\frac{n!}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}\prod\limits_{j=1}^{n}j^{m_{j}}}\left(\sum\limits_{m_{1}^{(1)}+\ldots+m_{1}^{(k)}=m_{1}}\frac{\theta_{1}^{m_{1}^{(1)}}\ldots\theta_{k}^{m_{1}^{(k)}}}{m_{1}^{(1)}!\ldots
m_{1}^{(k)}!}\right)\ldots\left(\sum\limits_{m_{n}^{(1)}+\ldots+m_{n}^{(k)}=m_{n}}\frac{\theta_{1}^{m_{n}^{(1)}}\ldots\theta_{k}^{m_{n}^{(k)}}}{m_{n}^{(1)}!\ldots
m_{n}^{(k)}!}\right)\\\
&=\frac{n!}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\frac{\left(\theta_{1}+\ldots+\theta_{k}\right)^{m_{1}+\ldots+m_{n}}}{\prod\limits_{j=1}^{n}j^{m_{j}}m_{j}!}.\end{split}$

The last expression is
$M_{\theta_{1}+\ldots+\theta_{k};n}^{\operatorname{Ewens}}\left(\lambda\right)$,
see equation (2.7). ∎

##  4\. The refined Ewens sampling formula as a distribution on a wreath
product. The Chinese restaurant process for wreath products

###  4.1. The refined Ewens sampling formula as a distribution on $G\sim S(n)$

Let $G$ be a finite group. Denote by $G_{\ast}$ the set of conjugacy classes
in $G$, and assume that $G_{\ast}$ consists of $k$ conjugacy classes, namely
$G_{\ast}=\left\\{c_{1},\ldots,c_{k}\right\\}$. Let $S(n)$ be the symmetric
group of degree $n$, i.e. the permutation group of
$\left\\{1,2,\ldots,n\right\\}$. The wreath product $G\sim S(n)$ is a group
whose underlying set is

$G^{n}\times
S(n)=\left\\{\left(\left(g_{1},\ldots,g_{n}\right),s\right):\;g_{i}\in G,s\in
S(n)\right\\}.$

The multiplication in $G\sim S(n)$ is defined by

$\left((g_{1},\ldots,g_{n}),s\right)\left((h_{1},\ldots,h_{n}),t\right)=\left((g_{1}h_{s^{-1}(1)},\ldots,g_{n}h_{s^{-1}(n)}),st\right).$

When $n=1$, $G\sim S(1)$ is $G$. The number of elements in $G\sim S(n)$ is
equal to $|G|^{n}n!$.

Let

(4.1) $x=\left(\left(g_{1},\ldots,g_{n}\right),s\right)\in G\sim S(n).$

The permutation $s$ can be written as a product of disjoint cycles. If
$(i_{1},i_{2},\ldots,i_{r})$ is one of these cycles, then the element
$g_{i_{r}}g_{i_{r-1}}\ldots g_{i_{1}}$ is called the cycle-product of
$(i_{1},i_{2},\ldots,i_{r})$. If $c_{l}$ is a conjugacy class in $G$, then a
cycle $(i_{1},\ldots,i_{r})$ of $s$ is called of type $c_{l}$ if its cycle-
product belongs to $c_{l}$. We denote by $[x]_{c_{l}}$ the number of cycles of
$s$ whose cycle-products belong to $c_{l}$.

Both the conjugacy classes and the irreducible representations of $G\sim S(n)$
are parameterized by multiple partitions
$\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$, where $k$
is the number of conjugacy classes in $G$, and where
$|\lambda^{(1)}|+\ldots+|\lambda^{(k)}|=n$, see, for example, Macdonald [15],
Appendix B. In particular, if
$x=\left(\left(g_{1},\ldots,g_{n}\right),s\right)$ belongs to the conjugacy
class of $G\sim S(n)$ parameterized by
$\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$, then

$\lambda^{(l)}=\left(1^{m_{1}^{(l)}(s)}2^{m_{2}^{(l)}(s)}\ldots
n^{m_{n}^{(l)}(s)}\right),$

where $m_{j}^{(l)}(s)$ is equal to the number of cycles of length $j$, and of
type $c_{l}$ in $s$.

To an element $x$ of $G\sim S(n)$ assign a multiple partition
$\Lambda_{n}^{(k)}$ labelling the conjugacy class of $x$. The projection

$x\longrightarrow\Lambda_{n}^{(k)}=\left(\lambda^{(1)},\ldots,\lambda^{(k)}\right)$

takes any central probability measure $P_{n}^{(k)}$ on $G\sim S(n)$ (central
measures are those invariant the conjugation action of $G\sim S(n)$) to a
probability measure $M_{n}^{(k)}$ on the set
$\operatorname{\mathbb{Y}}_{n}^{(k)}$ of all multiple partitions of $n$ into
$k$ components. This establishes a one-to-one correspondence

(4.2) $P_{n}^{(k)}\longleftrightarrow M_{n}^{(k)}$

between arbitrary central probability measures on $G\sim S(n)$, and arbitrary
probability measures on $\operatorname{\mathbb{Y}}_{n}^{(k)}$.

######  Proposition 4.1.

Under correspondence (4.2), the probability measure
$M^{\operatorname{Ewens}}_{\theta_{1},\ldots,\theta_{k}}$ defined by equation
(2.5) turns into a probability measure

(4.3)
$P_{t_{1},\ldots,t_{k};n}^{\operatorname{Ewens}}(x)=\frac{t_{1}^{[x]_{c_{1}}}t_{2}^{[x]_{c_{2}}}\ldots
t_{k}^{[x]_{c_{k}}}}{|G|^{n}\left(\frac{t_{1}}{\zeta_{c_{1}}}+\ldots+\frac{t_{k}}{\zeta_{c_{k}}}\right)_{n}}$

on $G\sim S(n)$. Here

$\zeta_{c_{l}}=\frac{|G|}{|c_{l}|},\;\;t_{l}=\frac{|G|}{|c_{l}|}\theta_{l},\;\;l=1,\ldots,k.$

###### Proof.

See Strahov [18], Section 4. ∎

######  Remark 4.2.

The probability measure $P_{t_{1},\ldots,t_{k};n}^{\operatorname{Ewens}}$ was
used in Strahov [19] to construct generalized regular representations of big
wreath products $G\sim S(\infty)$, and to develop harmonic analysis on $G\sim
S(\infty)$.

###  4.2. An analogue of the Chinese restaurant process

Our aim here is to construct a process that will generate elements of $G\sim
S(n)$ with the distribution defined by equation (4.3). This process will be an
analogue of the Chinese restaurant process [1], which is well-known in the
context of the symmetric group.

On the first step of the process one obtains an element $x=(g,(1))\in G\sim
S(1)$ with probability

$P^{\operatorname{Ewens}}_{t_{1},\ldots,t_{k};n=1}(x)=\left\\{\begin{array}[]{ll}\frac{t_{1}}{|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}},&g\in
c_{1},\\\ \vdots&\\\ \frac{t_{k}}{|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}},&g\in
c_{k}.\end{array}\right.$

Next, let us describe how an element $x^{\prime}\in G\sim S(n+1)$ is obtained
from a given element $x\in G\sim S(n)$ represented as in equation (4.1). Write
$s\in S(n)$ as a union of disjoint cycles,

(4.4) $s=\left(i_{1}i_{2}\ldots i_{m_{1}}\right)\left(j_{1}j_{2}\ldots
j_{m_{2}}\right)\ldots\left(\nu_{1}\nu_{2}\ldots\nu_{m_{p}}\right),$

and represent the elements $(g_{1},\ldots,g_{n})$ in (4.1) as a collection of
ordered lists (where each list corresponds to a relevant cycle of $s$)

(4.5)
$\left(g_{i_{1}},g_{i_{2}},\ldots,g_{i_{m_{1}}}\right)\left(g_{j_{1}},g_{j_{2}},\ldots,g_{j_{m_{2}}}\right)\ldots\left(g_{\nu_{1}},g_{\nu_{2}},\ldots,g_{\nu_{m_{p}}}\right).$

Clearly, $x\in G\sim S(n)$ represented as in (4.1) determines (4.4), (4.5)
uniquely. Conversely, collection (4.4) of cycles, and collection (4.5) of
ordered lists of elements of $G$ gives us unambiguously an element $x$ of
$G\sim S(n)$.

The process is such that an element $x\in G\sim S(n)$ can be lifted to an
element of $G\sim S(n+1)$ in two ways. First, we can insert $n+1$ to an
existing cycle of $s$, and change the corresponding lists of elements of $G$
accordingly. If

$\left(i_{1}i_{2}\ldots i_{s}i_{s+1}\ldots
i_{m_{1}}\right)\rightarrow\left(i_{1}i_{2}\ldots i_{s}\;n+1\;i_{s+1}\ldots
i_{m_{1}}\right),$

then we replace

(4.6) $\left(g_{i_{1}},g_{i_{2}},\ldots,g_{i_{s}},g_{i_{s+1}},\ldots
g_{i_{m_{1}}}\right)$

by

(4.7)
$\left(g_{i_{1}},g_{i_{2}},\ldots,g_{i_{s}}\left(g_{n+1}\right)^{-1},g_{n+1},g_{i_{s+1}},\ldots
g_{i_{m_{1}}}\right).$

In this way a new element $x^{\prime}\in G\sim S(n+1)$ will be obtained. Note
that the condition $[x^{\prime}]_{c_{l}}=[x]_{c_{l}}$ is satisfied for
$l=1,\ldots,k$ in this case. The probability that $x^{\prime}\in G\sim S(n+1)$
is obtained from a given $x\in G\sim S(n)$ as described above is assigned to
be equal to

$\frac{1}{|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}+n|G|}.$

Second, and integer $n+1$ can start a new cycle. In this case the permutation
$s$ is replaced by

(4.8) $s^{\prime}=\left(i_{1}i_{2}\ldots
i_{m_{1}}\right)\left(j_{1}j_{2}\ldots
j_{m_{2}}\right)\ldots\left(\nu_{1}\nu_{2}\ldots\nu_{m_{p}}\right)(n+1)\in
S(n+1).$

In words, $s^{\prime}\in S(n+1)$ is obtained from $s\in S(n)$ by attaching the
new cycle $(n+1)$. The collection of ordered lists (4.5) is replaced by

(4.9)
$\left(g_{i_{1}},g_{i_{2}},\ldots,g_{i_{m_{1}}}\right)\left(g_{j_{1}},g_{j_{2}},\ldots,g_{j_{m_{2}}}\right)\ldots\left(g_{\nu_{1}},g_{\nu_{2}},\ldots,g_{\nu_{m_{p}}}\right)(g_{n+1}).$

Expressions (4.8) and (4.9) determine an element $x^{\prime}\in G\sim S(n+1)$.
If $g_{n+1}\in c_{l}$, then $[x^{\prime}]_{c_{l}}=[x]_{c_{l}}+1$, and the
probability of the event that $x^{\prime}$ is obtained from $x$ is assigned to
be equal to

$\frac{t_{l}}{|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}+n|G|}.$

If $G$ contains the unit element only, then the wreath product $G\sim S(n)$ is
isomorphic to the symmetric group $S(n)$, and the process described above
coincides with the Chinese restaurant process. For this reason we will refer
to the process described above as to the generalized Chinese restaurant
process.

######  Proposition 4.3.

The generalized Chinese restaurant process generates the distribution on
$G\sim S(n)$ given by formula (4.3).

###### Proof.

The probability measure $P_{t_{1},\ldots,t_{k};n}^{\operatorname{Ewens}}$ on
$G\sim S(n)$ defined by equation (4.3) can be represented as

(4.10) $\begin{split}&P_{t_{1},\ldots,t_{k};n}^{\operatorname{Ewens}}(x)\\\
&=\frac{t_{1}^{[x]_{c_{1}}}t_{2}^{[x]_{c_{2}}}\ldots
t_{k}^{[x]_{c_{k}}}}{\left(|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}\right)\left(|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}+|G|\right)\ldots\left(|c_{1}|t_{1}+\ldots+|c_{k}|t_{k}+(n-1)|G|\right)}.\end{split}$

It is straightforward to check that the generalized Chinese restaurant process
described above generates the distribution on $G\sim S(n)$ given by formula
(4.10). ∎

##  5\. Numbers of alleles

Let $K_{n}^{(l)}$, $l\in\\{1,\ldots,k\\}$, be the number of different alleles
of class $l$ in a random sample of size $n$. The purpose of this Section is to
derive exact formulae for the random variables $K_{n}^{(l)}$, and to describe
their limiting distributions in different asymptotic regimes under the
assumption that the allelic composition of the sample is governed by the
refined Ewens sampling formula, equation (1.3).

###  5.1. Exact formulas for $K_{n}^{(l)}$

In the context of alleles, the sampling of a black object of the $l$th class
in the generalized Hoppe urn described in Section 3.2 is a mutation
characterized by the parameter $\theta_{l}$ (where $l\in\\{1,\ldots,k\\}$),
which induces new allelic types of class $l$. The selection of a non-black
object of the $l$th class has the interpretation that the genes are
multiplying themselves within their allelic type. With this interpretation, we
have

(5.1) $K_{n}^{(l)}=Y_{1}^{(l)}+Y_{2}^{(l)}+\ldots+Y_{n}^{(l)},$

where $Y_{j}^{(l)}$ is the indicator that assumes $1$ if the black object of
class $l$ (with mass $\theta_{l}>0$) is chosen in the $j$th draw, and assumes
$0$ otherwise. Then we can write

(5.2)
$\operatorname{Prob}\left(Y_{j}^{(l)}=1\right)=\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+j-1},\;\;\operatorname{Prob}\left(\operatorname{\mathbb{Y}}_{j}^{(l)}=0\right)=1-\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+j-1}.$

Note that the Bernoulli random variables
$\operatorname{\mathbb{Y}}_{1}^{(l)}$,
$\ldots$,$\operatorname{\mathbb{Y}}_{n}^{(l)}$ are independent. Taking
expectations we find

(5.3)
$\mathbb{E}\left(K_{n}^{(l)}\right)=\mathbb{E}(Y_{1}^{(l)})+\mathbb{E}(Y_{2}^{(l)})+\ldots+\mathbb{E}(Y_{n}^{(l)}),$

which gives

(5.4)
$\mathbb{E}\left(K_{n}^{(l)}\right)=\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}}+\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+1}+\ldots+\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+n-1}.$

Denote by $\mathcal{H}_{n}^{(p)}(x)$ a function defined by

(5.5)
$\mathcal{H}_{n}^{(p)}(x)=\frac{1}{x^{p}}+\frac{1}{(x+1)^{p}}+\ldots+\frac{1}{(x+n-1)^{p}},$

where $p\in\\{1,2,\ldots\\}$. Then we obtain

(5.6)
$\mathbb{E}\left(K_{n}^{(l)}\right)=\theta_{l}\mathcal{H}_{n}^{(1)}\left(\theta_{1}+\ldots+\theta_{k}\right).$

Independence of the Bernoulli random variables
$\operatorname{\mathbb{Y}}_{1}^{(l)}$,
$\ldots$,$\operatorname{\mathbb{Y}}_{n}^{(l)}$ implies

(5.7)
$\operatorname{Var}\left(K_{n}^{(l)}\right)=\operatorname{Var}(Y_{1}^{(l)})+\operatorname{Var}(Y_{2}^{(l)})+\ldots+\operatorname{Var}(Y_{n}^{(l)}),$

which gives

(5.8)
$\operatorname{Var}\left(K_{n}^{(l)}\right)=\sum\limits_{j=1}^{n}\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+j-1}\left(1-\frac{\theta_{l}}{\theta_{1}+\ldots+\theta_{k}+j-1}\right).$

Therefore,

(5.9)
$\operatorname{Var}\left(K_{n}^{(l)}\right)=\theta_{l}\mathcal{H}_{n}^{1)}\left(\theta_{1}+\ldots+\theta_{k}\right)-\theta_{l}^{2}\mathcal{H}_{n}^{(2)}\left(\theta_{1}+\ldots+\theta_{k}\right).$

Denote by $\left[\begin{array}[]{c}n\\\ m\end{array}\right]$ be the number of
permutations from $S(n)$ with $m$ disjoint cycles. This number can be written
explicitly as

$\left[\begin{array}[]{c}n\\\
m\end{array}\right]=\underset{l(\lambda)=m}{\sum\limits_{\lambda\vdash
n}}\frac{n!}{\prod_{j=1}^{n}j^{m_{j}(\lambda)}\left(m_{j}(\lambda)\right)!}.$

With this notation we can write the joint distribution of the random variables
$K_{n}^{(1)}$, $\ldots$, $K_{n}^{(k)}$ as

(5.10)
$\operatorname{Prob}\left(K_{n}^{(1)}=p_{1},\ldots,K_{n}^{(k)}=p_{k}\right)=\frac{n!\theta_{1}^{p_{1}}\ldots\theta_{k}^{p_{k}}}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\sum\limits_{n_{1}+\ldots+n_{k}=n}\frac{\left[\begin{array}[]{c}n_{1}\\\
p_{1}\end{array}\right]\ldots\left[\begin{array}[]{c}n_{k}\\\
p_{k}\end{array}\right]}{n_{1}!\ldots n_{k}!}.$

In the case of $k=1$ (all alleles belong to the same class) formula (5.4) for
the expectation of $K_{n}^{(l)}$, formula (5.9) for the variance of
$K_{n}^{(l)}$, and formula (5.10) for the joint distribution of random
variables $K_{n}^{(1)}$, $\ldots$, $K_{n}^{(k)}$ are reduced to the well-known
formulae derived by Ewens in Ref. [7].

###  5.2. Inequalities for $\mathcal{H}_{n}^{(1)}(x)$ and
$\mathcal{H}_{n}^{(2)}(x)$

In the subsequent asymptotic analysis we will need the following inequalities
for the functions $\mathcal{H}_{n}^{(1)}(x)$ and $\mathcal{H}_{n}^{(2)}(x)$
defined by equation (5.5).

######  Lemma 5.1.

Let the functions $\mathcal{H}_{n}^{(1)}(x)$ and $\mathcal{H}_{n}^{(2)}(x)$ be
defined by equation (5.5). Then the following inequalities hold true

(5.11)
$\frac{n}{2x(x+n)}<\mathcal{H}_{n}^{(1)}(x)-\log\left(1+\frac{n}{x}\right)<\frac{n}{x(x+n)},$

and

(5.12)
$\frac{n}{x(x+n)}<\mathcal{H}_{n}^{(2)}(x)<\frac{1}{x^{2}}+\frac{n-1}{x(x+n-1)}.$

Here $x>0$, and $n=1,2,\ldots$.

###### Proof.

The proof of these inequalities can be extracted from that of Proposition 1 in
Tsukuda [21]. We reproduce the proof here for reader’s convenience.

In order to prove (5.1) note that

$\mathcal{H}_{n}^{(1)}(x)-\log\left(1+\frac{n}{x}\right)=\sum\limits_{j=1}^{n}\left(\frac{1}{x+j-1}-\int\limits_{x+j-1}^{x+j}\frac{dt}{t}\right).$

Since

$\frac{1}{x+j}<\int\limits_{x+j-1}^{x+j}\frac{dt}{t}<\frac{1}{x+j-1},\;\;\;\;\int\limits_{x+j-1}^{x+j}\frac{dt}{t}<\frac{1}{2}\left(\frac{1}{x+j-1}+\frac{1}{x+j}\right),$

we have

$\frac{1}{2}\left(\frac{1}{x+j-1}-\frac{1}{x+j}\right)<\frac{1}{x+j-1}-\int\limits_{x+j-1}^{x+j}\frac{dt}{t}<\frac{1}{x+j-1}-\frac{1}{x+j}.$

This gives

$\frac{1}{2}\sum\limits_{j=1}^{n}\left(\frac{1}{x+j-1}-\frac{1}{x+j}\right)<\mathcal{H}_{n}^{(1)}(x)-\log\left(1+\frac{n}{x}\right)<\sum\limits_{j=1}^{n}\left(\frac{1}{x+j-1}-\frac{1}{x+j}\right),$

and inequality (5.1) follows. In order to prove (5.12) note that

$\sum\limits_{j=1}^{n}\frac{1}{\left(x+j-1\right)^{2}}>\int\limits_{x}^{x+1}\frac{dt}{t^{2}}=\frac{n}{x(x+n)},$

and that

$\frac{n-1}{x(x+n-1)}=\int\limits_{x}^{x+n-1}\frac{dt}{t^{2}}>\frac{1}{(x+1)^{2}}+\ldots+\frac{1}{(x+n-1)^{2}}=\sum\limits_{j=1}^{n}\frac{1}{(x+j-1)^{2}}-\frac{1}{x^{2}}.$

∎

###  5.3. Limit theorems

####  5.3.1. The case of constant mutation parameters

First, let us assume that the mutation parameters $\theta_{1}$, $\ldots$,
$\theta_{k}$ are constant strictly positive numbers independent on the sample
size $n$. Then the following theorem can be easily obtained from the results
of Sections 5.1 and 5.2

######  Theorem 5.2.

Let $K_{n}^{(l)}$, $l\in\\{1,\ldots,k\\}$, be the number of different alleles
of class $l$ in a random sample of size $n$, and assume that the mutation
rates $\theta_{1}$, $\ldots$, $\theta_{k}$ are constant strictly positive
numbers. Then for each $l$, $l\in\left\\{1,\ldots,k\right\\}$, we have

(5.13) $\mathbb{E}\left(K_{n}^{(l)}\right)=\theta_{l}\log(n)+O(1),$ (5.14)
$\operatorname{Var}\left(K_{n}^{(l)}\right)=\theta_{l}\log(n)+O(1),$ (5.15)
$\frac{K_{n}^{(l)}}{\log(n)}\overset{P}{\longrightarrow}\theta_{l},$

and

(5.16)
$\frac{K_{n}^{(l)}-\theta_{l}\log(n)}{\sqrt{\log(n)}}\overset{D}{\longrightarrow}\mathcal{N}(0,1),$

as $n\rightarrow\infty$. In (5.15) the convergence is in probability, in
(5.16) the convergence is in distribution. We have denoted by
$\mathcal{N}(0,1)$ a standard normal random variable.

###### Proof.

Formulae (5.13) and (5.14) follow from (5.6), (5.7), (5.1) and (5.12). The
proof of (5.15) and (5.16) is a repetition of arguments for the standard Ewens
distribution, see the book by Mahmoud [16], Section 9. ∎

######  Remark 5.3.

The asymptotic normality of $K_{n}^{(l)}$ (equation (5.16)) is an extension of
the Watterson limit theorem [23] to the case of an arbitrary number $k$ of
classes of alleles.

####  5.3.2. Growth of mutation parameters with $n$

It is of interest to investigate asymptotic properties of numbers of
$K_{n}^{(l)}$ in the situation where the mutation parameters $\theta_{1}$,
$\ldots$, $\theta_{k}$ tends to infinity as the size $n$ of a sample
increases. In the case of $k=1$ (i.e.e in the case of the classical Ewens
formula) different asymptotic regimes where both $n$ and the mutation
parameter tend to infinity were considered in Feng [11], Tsukuda [21, 22]. Our
aim here is to extend some results obtained in Ref. [21] to the situation
where the allelic composition of a sample is described by the refined Ewens
sampling formula, equation (1.3).

######  Proposition 5.4.

Let $K_{n}^{(l)}$, $l\in\\{1,\ldots,k\\}$, be the number of different alleles
of class $l$ in a random sample of size $n$. If the mutation parameters
$\theta_{1}$, $\ldots$, $\theta_{k}$ are strictly positive valued functions of
$n$ growing to infinity as $n\rightarrow\infty$, then

(5.17)
$\mathbb{E}\left(K_{n}^{(l)}\right)=\theta_{l}\log\left(1+\frac{n}{w}\right)+O(1),$

and

(5.18)
$\operatorname{Var}\left(K_{n}^{(l)}\right)=\theta_{l}\log\left(1+\frac{n}{w}\right)-\frac{n\theta_{l}^{2}}{w(w+n)}+O(1),$

as $n\rightarrow\infty$, where $w$ is the sum of the mutation parameters,

(5.19) $w=\theta_{1}+\ldots+\theta_{k}.$

###### Proof.

Formula (5.6) and inequality (5.11) imply

(5.20)
$\frac{\theta_{l}n}{2w(w+n)}<\mathbb{E}\left(K_{n}^{(l)}\right)-\theta_{l}\log\left(1+\frac{n}{w}\right)<\frac{\theta_{l}n}{w(w+n)},$

which gives equation (5.17). From inequality (5.12) we conclude that

(5.21)
$0<\theta_{l}^{2}\mathcal{H}_{n}^{(2)}(w)-\frac{n\theta^{2}_{l}}{w(w+n)}<\frac{\theta_{l}^{2}}{w^{2}}+\frac{\theta_{l}^{2}}{(w+n)(w+n-1)}.$

Therefore, equation (5.18) follows from (5.6),(5.9), (5.17), and (5.21). ∎

Next, we investigate the convergence in probability of random variables
$K_{n}^{(l)}$, $l\in\left\\{1,\ldots,k\right\\}$.

######  Theorem 5.5.

Assume that the mutation parameters $\theta_{1}$, $\ldots$, $\theta_{k}$
depend on $n$ as

(5.22)
$\theta_{1}(n)=\alpha_{1}n^{\beta},\ldots,\theta_{k}(n)=\alpha_{k}n^{\beta};\;\;\beta\geq
0;\;\;\alpha_{1}>0,\ldots,\alpha_{k}>0.$

Let $K_{n}^{(l)}$, $l\in\\{1,\ldots,k\\}$, be the number of different alleles
of class $l$ in a random sample of size $n$.
(a) If $0\leq\beta<1$, then

$\frac{K_{n}^{(l)}}{n^{\beta}\log(n)}\overset{P}{\longrightarrow}\alpha_{l}(1-\beta).$

(b) If $\beta=1$, then

$\frac{K_{n}^{(l)}}{n}\overset{P}{\longrightarrow}\alpha_{l}\log\frac{1+A}{A},$

where

(5.23) $A=\alpha_{1}+\ldots+\alpha_{k}.$

(c) If $\beta>1$, then

$\frac{K_{n}^{(l)}}{n}\overset{P}{\longrightarrow}\frac{\alpha_{l}}{A},$

where $A$ is defined by equation (5.23).

###### Proof.

If the dependence of the mutation parameters on $n$ is giving by equation
(5.22), then the following asymptotic expressions for
$\mathbb{E}\left(K_{n}^{(l)}\right)$ and
$\operatorname{Var}\left(K_{n}^{(l)}\right)$ can be obtained from equations
(5.17) and (5.18).
(a) If $0\leq\beta<1$, then

$\mathbb{E}\left(K_{n}^{(l)}\right)=\alpha_{l}(1-\beta)n^{\beta}\log(n)+O\left(n^{\beta}\right),$

and

$\operatorname{Var}\left(K_{n}^{(l)}\right)=\alpha_{l}(1-\beta)n^{\beta}\log(n)+O\left(n^{\beta}\right),$

as $n\rightarrow\infty$.
(b) If $\beta=1$, then

$\mathbb{E}\left(K_{n}^{(l)}\right)=\alpha_{l}n\log\left(\frac{1+A}{A}\right)+O\left(1\right),$

where $A$ is defined by equation (5.23), and

$\operatorname{Var}\left(K_{n}^{(l)}\right)=\alpha_{l}n\log\left(\frac{1+A}{A}\right)-\frac{\alpha_{l}^{2}n}{A(1+A)}+O\left(1\right),$

as $n\rightarrow\infty$.
(c) If $\beta>1$, then

$\mathbb{E}\left(K_{n}^{(l)}\right)=\frac{\alpha_{l}}{A}n+O\left(n^{2-\beta}\right),$

and

$\operatorname{Var}\left(K_{n}^{(l)}\right)=\frac{\alpha_{l}}{A}\left(1-\frac{\alpha_{l}}{A}\right)n+O\left(n^{2-\beta}\right),$

as $n\rightarrow\infty$. Note that for each $\beta\geq 0$ the condition

(5.24)
$\frac{\operatorname{Var}\left(K_{n}^{(l)}\right)}{\left(\mathbb{E}\left(K_{n}^{(l)}\right)\right)^{2}}\longrightarrow
0,\;\;\mbox{as}\;\;n\longrightarrow\infty,$

is satisfied. Recall that by Chebyshev’s inequality, for any fixed
$\epsilon>0$

$\operatorname{Prob}\left(\left|K_{n}^{(l)}-\mathbb{E}\left(K_{n}^{(l)}\right)\right|>\epsilon\right)\leq\frac{\operatorname{Var}\left(K_{n}^{(l)}\right)}{\epsilon^{2}}.$

Replacing $\epsilon$ by $\epsilon\mathbb{E}\left(K_{n}^{(l)}\right)$, we
obtain

(5.25)
$\operatorname{Prob}\left(\left|\frac{K_{n}^{(l)}}{\mathbb{E}\left(K_{n}^{(l)}\right)}-1\right|>\epsilon\right)\leq\frac{\operatorname{Var}\left(K_{n}^{(l)}\right)}{\epsilon^{2}\left(\mathbb{E}\left(K_{n}^{(l)}\right)\right)^{2}}.$

It follows from (5.24) and (5.25) that

(5.26)
$\frac{K_{n}^{(l)}}{\mathbb{E}\left(K_{n}^{(l)}\right)}\overset{P}{\longrightarrow}1,$

as $n\rightarrow\infty$, for each $\beta\geq 0$. The results in the statement
of the Theorem follow from (5.26), and from the convergence relations

(5.27)
$\frac{\mathbb{E}\left(K_{n}^{(l)}\right)}{n^{\beta}\log(n)}\longrightarrow\alpha_{l}(1-\beta),\;\;\;0\leq\beta<1,$
(5.28)
$\frac{\mathbb{E}\left(K_{n}^{(l)}\right)}{n}\longrightarrow\alpha_{l}\log\left(\frac{1+A}{A}\right)-\frac{\alpha_{l}^{2}}{A(1+A)},\;\;\;\beta=1,$

and

(5.29)
$\frac{\mathbb{E}\left(K_{n}^{(l)}\right)}{n}\longrightarrow\frac{\alpha_{l}}{A},\;\;\;\beta>1,$

as $n\rightarrow\infty$. ∎

######  Remark 5.6.

Equation (5.24) leads to the fact that $K_{n}^{(l)}$ is highly concentrated
around its expectation in the case were the mutation parameters depend on $n$
as in equation (5.22). However, if the mutation parameters are allowed to grow
with different rates, then (5.24) might not be true. For example, if $k=2$,
$\theta_{1}(n)=\alpha_{1}n^{1/2}$, and $\theta_{2}(n)=\alpha_{2}n^{3/2}$, then
we obtain that $\mathbb{E}\left(K_{n}^{(l)}\right)=O(1)$,
$\operatorname{Var}\left(K_{n}^{(l)}\right)=O(1)$, and (5.24) is not valid.

If $k=1$, the asymptotic normality of $K_{n}^{(l)}$ as $n\rightarrow\infty$
(with a dependence of the mutation parameter on $n$) was established in
Tsukuda [21]. Theorem 5.7 below presents the asymptotic normality of
$K_{n}^{(l)}$ when the number $k$ of classes of alleles is arbitrary, and the
mutation parameters $\theta_{1}$, $\ldots$, $\theta_{k}$ depend on a sample
size $n$ as in equation (5.22).

######  Theorem 5.7.

Let $K_{n}^{(l)}$, $l\in\\{1,\ldots,k\\}$, be the number of different alleles
of class $l$ in a random sample of size $n$, and assume that the mutation
parameters are given by equation (5.22).
(a) If $0<\beta\leq\frac{3}{2}$, then for each
$l\in\left\\{1,\ldots,k\right\\}$

(5.30)
$\frac{K_{n}^{(l)}-\theta_{l}\log\left(1+\frac{n}{w}\right)}{\sqrt{\theta_{l}\log\left(1+\frac{n}{w}\right)-\frac{n\theta_{l}^{2}}{w\left(w+n\right)}}}\overset{D}{\longrightarrow}\mathcal{N}\left(0,1\right),$

as $n\rightarrow\infty$, where $w$ is the sum of the mutation parameters
defined by equation (5.19), and $\mathcal{N}(0,1)$ denotes a standard normal
random variable.
(b) If $\beta>\frac{3}{2}$ and $k\geq 2$, then

(5.31)
$\frac{K_{n}^{(l)}-n\frac{\theta_{l}}{w}}{\sqrt{n\frac{\theta_{l}}{w}\left(1-\frac{\theta_{l}}{w}\right)}}\overset{D}{\longrightarrow}\mathcal{N}(0,1)$

is satisfied for each $l\in\left\\{1,\ldots,k\right\\}$, as
$n\rightarrow\infty$.

######  Remark 5.8.

(a) If $k=1$, then $\frac{\theta_{l}}{w}=1$, and (5.31) does not hold true.
For the results on $K_{n}^{(l)}$ in the case $k=1$ we refer the reader to
Theorem 2 and Theorem 3 in Tsukuda [21].
(b) If $\beta=0$, then the convergence in distribution to a standard normal
random variable is given by equation (5.16).
(c) If $k=1$, and $n^{2}/\theta_{1}\rightarrow c>0$ as $n\rightarrow\infty$,
then Tsukuda [21] showed that

(5.32)
$\frac{K_{n}^{(l)}-\theta_{1}\log\left(1+\frac{n}{\theta_{1}}\right)}{\sqrt{\theta_{1}\left(\log\left(1+\frac{n}{\theta_{1}}\right)-\frac{n}{n+\theta_{1}}\right)}}\overset{D}{\longrightarrow}\frac{c/2-P_{c/2}}{\sqrt{c/2}},$

where $P_{c/2}$ is a Poisson variable with mean $c/2$. However, Theorem 5.7,
(b) says that for $k\geq 2$ we still have convergence in distribution to a
standard normal random variable.

###### Proof.

The proof here is an extension of that in Tsukuda [21], Section 3.2 to the
case of an arbitrary number $k$ of mutation parameters.
(a) Assume that $0<\beta\leq\frac{3}{2}$. Set

(5.33)
$\Delta_{l}=\theta_{l}\log\left(1+\frac{n}{w}\right),\;\;\sigma_{l}^{2}=\theta_{l}\log\left(1+\frac{n}{w}\right)-\frac{n\theta_{l}^{2}}{w(w+n)}.$

We note that

(5.34) $\Delta_{l}=\sigma_{l}^{2}+\frac{n\theta_{l}^{2}}{w(w+n)}.$

Using (5.1) and independence of the Bernoulli random variables $Y_{1}^{(l)}$,
$\ldots$, $Y_{n}^{(l)}$ we find

(5.35)
$\mathbb{E}\left(e^{\frac{K_{n}^{(l)}-\Delta_{l}}{\sigma_{l}}t}\right)=e^{-\frac{\Delta_{l}}{\sigma_{l}}t}\frac{\Gamma\left(\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w+n\right)}{\Gamma\left(w+n\right)}\frac{\Gamma(w)}{\Gamma\left(\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w\right)}.$

Application of the Stirling formula gives

(5.36)
$\begin{split}&\frac{\Gamma\left(\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w+n\right)}{\Gamma\left(w+n\right)}\frac{\Gamma(w)}{\Gamma\left(\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w\right)}\\\
&\sim\frac{\left(\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w+n\right)^{\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w+n-1}}{(w+n)^{w+n-1}}\;\frac{w^{w-1}}{\left(\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w\right)^{\theta_{l}\left(e^{\frac{t}{\sigma_{l}}}-1\right)+w-1}},\end{split}$

as $n\rightarrow\infty$. The right-hand side of (5.36) can be rewritten as

(5.37)
$\begin{split}&\left[1+\frac{n}{w}\right]^{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}\right]^{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+w+n-1}\\\
&\times\left[1+\frac{\theta_{l}}{w}\left(e^{t/\sigma_{l}}-1\right)\right]^{-\theta_{l}\left(e^{t/\sigma_{l}}-1\right)-w+1}.\end{split}$

The equations (5.35), (5.36), (5.33)-(5.37) imply

(5.38)
$\begin{split}&\log\left[\mathbb{E}\left(e^{\frac{K_{n}^{(l)}-\Delta_{l}}{\sigma_{l}}t}\right)\right]\sim-\sigma_{l}t-\frac{n\theta_{l}^{2}t}{\sigma_{l}w(w+n)}+\sigma_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)\\\
&+\frac{n\theta_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)}{w(w+n)}+\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+n-1+w\right]\log\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}\right]\\\
&-\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+w-1\right]\log\left[1+\frac{\theta_{l}}{w}\left(e^{t/\sigma_{l}}-1\right)\right].\\\
&=-\sigma_{l}t-\frac{n\theta_{l}^{2}t}{\sigma_{l}w(w+n)}+\sigma_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)\\\
&+\frac{n\theta_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)}{w(w+n)}+\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+n+w\right]\log\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}\right]\\\
&-\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+w\right]\log\left[1+\frac{\theta_{l}}{w}\left(e^{t/\sigma_{l}}-1\right)\right]+o(1),\end{split}$

as $n\rightarrow\infty$.

Let us consider the case where $0<\beta<3/2$. Below it is important that if
$0<\beta<3/2$, then $\theta_{l}/\sigma_{l}^{3}\rightarrow 0$ as
$n\rightarrow\infty$. This follows from equation (5.33), from (5.18), and from
explicit asymptotic formulae for the variance of $K_{n}^{(l)}$ in the proof of
Theorem 5.5. We can write

(5.39)
$\begin{split}&\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+n-1+w\right]\log\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}\right]\\\
&=\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+n+w\right]\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}-\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+n+w\right]\frac{\theta_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)^{2}}{2(w+n)^{2}}+o(1)\\\
&=\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+\frac{1}{2}\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)^{2}}{w+n}+o(1),\end{split}$

since $\frac{\theta_{l}^{3}(e^{t/\sigma_{l}}-1)^{3}}{2(w+n)^{2}}=o(1)$
provided that $\theta_{l}/\sigma_{l}^{3}\rightarrow 0$ as
$n\rightarrow\infty$. Also,

(5.40)
$\begin{split}&-\left[\theta_{l}\left(e^{t/\sigma_{l}}-1\right)+w-1\right]\log\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w}\right]\\\
&=-\theta_{l}\left(e^{t/\sigma_{l}}-1\right)-\frac{1}{2}\frac{\theta_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)^{2}}{w}+o(1),\end{split}$

as $n\rightarrow\infty$. Taking (5.39), (5.40) into account we rewrite (5.38)
as

(5.41)
$\begin{split}&\log\left[\mathbb{E}\left(e^{\frac{K_{n}^{(l)}-\Delta_{l}}{\sigma_{l}}t}\right)\right]\sim-\sigma_{l}t+\sigma_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)-\frac{n\theta_{l}^{2}t}{\sigma_{l}w(w+n)}\\\
&+\frac{n\theta^{2}_{l}\left(e^{t/\sigma_{l}}-1\right)}{w(w+n)}-\frac{n\theta^{2}_{l}\left(e^{t/\sigma_{l}}-1\right)^{2}}{2w(w+n)}.\end{split}$

The first two terms in the right-hand side of equation (5.41) give

(5.42)
$-\sigma_{l}t+\sigma_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)=\frac{t^{2}}{2}+o(1),$

as $n\rightarrow\infty$. It is easy to check that the last three terms in the
right-hand side of equation (5.41) are equal to $o(1)$, as
$n\rightarrow\infty$. We conclude that

(5.43)
$\log\left[\mathbb{E}\left(e^{\frac{K_{n}^{(l)}-\Delta_{l}}{\sigma_{l}}t}\right)\right]\sim\frac{t^{2}}{2},$

as $n\rightarrow\infty$. This implies the statement of Theorem 5.7 in the case
$0<\beta<3/2$.

If $\beta=3/2$, then we rewrite the two last terms in the right-hand side of
(5.38) as

(5.44)
$\begin{split}&(w+n)\log\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}\right]-w\log\left[1+\frac{\theta_{l}}{w}\left(e^{t/\sigma_{l}}-1\right)\right]\\\
&+\theta_{l}\left(e^{t/\sigma_{l}}-1\right)\log\left[1+\frac{\theta_{l}\left(e^{t/\sigma_{l}}-1\right)}{w+n}\right]-\theta_{l}\left(e^{t/\sigma_{l}}-1\right)\log\left[1+\frac{\theta_{l}}{w}\left(e^{t/\sigma_{l}}-1\right)\right]\\\
\end{split}$

As $n\rightarrow\infty$, the first two terms in the expression above are equal
to

$\frac{n\theta_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)^{2}}{2w(w+n)}+o(1),$

and the third and the fourth terms in (5.44) are equal to

$-\frac{n\theta_{l}^{2}\left(e^{t/\sigma_{l}}-1\right)^{2}}{w(w+n)}+o(1).$

Therefore, equations (5.41)-(5.43) are still hold true in the case
$\beta=3/2$, which completes the proof of Theorem 5.7, (a).

(b) Assume that $\beta>3/2$. We use the formula

(5.45)
$\begin{split}&\mathbb{E}\left(e^{\frac{K_{n}^{(l)}-\tilde{\Delta_{l}}}{\tilde{\sigma}_{l}}t}\right)\sim
e^{-\frac{\tilde{\Delta_{l}}}{\tilde{\sigma}_{l}}t}\frac{\left(\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w+n\right)^{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w+n-1}}{(w+n)^{w+n-1}}\\\
&\times\frac{w^{w-1}}{\left(\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w\right)^{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w-1}},\end{split}$

where

(5.46)
$\tilde{\Delta_{l}}=\frac{\theta_{l}}{w},\;\;\;\tilde{\sigma}_{l}^{2}=n\frac{\theta_{l}}{w}\left(1-\frac{\theta_{l}}{w}\right).$

We note that

(5.47)
$\begin{split}&\frac{\left[\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w+n\right]^{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w+n-1}}{\left[\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w\right]^{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w-1}}\\\
&=\left[1+\frac{n}{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w}\right]^{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w-1}\left[\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w+n\right]^{n}.\end{split}$

Therefore, we can rewrite equation (5.45) as

(5.48)
$\begin{split}&\mathbb{E}\left(e^{\frac{K_{n}^{(l)}-\tilde{\Delta_{l}}}{\tilde{\sigma}_{l}}t}\right)\sim
e^{-\frac{\tilde{\Delta_{l}}}{\tilde{\sigma}_{l}}t}\left[1+\frac{n}{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w}\right]^{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w-1}\\\
&\times\left[1+\frac{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)}{w+n}\right]^{n}\left(1+\frac{n}{w}\right)^{w-1}.\end{split}$

The crucial observation is that $\tilde{\sigma}_{l}\sim n^{1/2}$, and that

(5.49)
$n\log\left[1+\frac{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)}{w+n}\right]=n\frac{\theta_{l}t}{w\tilde{\sigma}_{l}}+n\frac{\theta_{l}t^{2}}{2w\tilde{\sigma}_{l}^{2}}-\frac{n\theta_{l}^{2}t^{2}}{2w^{2}\tilde{\sigma}_{l}^{2}}+o(1)=\frac{\tilde{\Delta_{l}}t}{\tilde{\sigma}_{l}}+\frac{t^{2}}{2}+o(1),$

as $n\rightarrow\infty$. In addition, it is easy to see that under the
assumption $\beta>3/2$

(5.50)
$\left[\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w-1\right]\log\left[1+\frac{n}{\theta_{l}\left(e^{t/\tilde{\sigma}_{l}}-1\right)+w}\right]-(w-1)\log\left[1+\frac{n}{w}\right]=o(1),$

as $n\rightarrow\infty$. Therefore, (5.43) holds true, which implies (5.31).
Theorem 5.7 is proved. ∎

##  6\. Representation in terms of random matrices with Poisson entries

Denote by $A(n)$ a random matrix of size $n\times k$ whose entries
$A_{j}^{(l)}(n)$ (where $j=1,\ldots,n$, and $l=1,\ldots,k$) are random
variables taking values in the set $\left\\{0,1,2,\ldots\right\\}$ of positive
integers. Assume that the joint probability distribution of $A_{j}^{(l)}(n)$
is given by the refined Ewens sampling formula (equation (1.3)). Then
$A_{j}^{(l)}(n)$ represents the number of alleles of class $l$ that appear $j$
times in the sample of size $n$.

######  Proposition 6.1.

Let

(6.1) $\left(\begin{array}[]{ccc}a_{1}^{(1)}&\ldots&a_{1}^{(k)}\\\ \vdots&&\\\
a_{n}^{(1)}&\ldots&a_{n}^{(k)}\end{array}\right),$

be a fixed matrix whose entries are non-negative integers satisfying the
condition

(6.2) $\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}ja_{j}^{(l)}=n.$

Then for each $n\geq 2$

(6.3)
$\begin{split}&\operatorname{Prob}\left(A(n)=\left(\begin{array}[]{ccc}a_{1}^{(1)}&\ldots&a_{1}^{(k)}\\\
\vdots&&\\\ a_{n}^{(1)}&\ldots&a_{n}^{(k)}\end{array}\right)\right)\\\
&=\operatorname{Prob}\left(\left(\begin{array}[]{ccc}\eta_{1}(\theta_{1})&\ldots&\eta_{1}(\theta_{k})\\\
\vdots&&\\\
\eta_{n}(\theta_{1})&\ldots&\eta_{n}(\theta_{k})\end{array}\right)=\left(\begin{array}[]{ccc}a_{1}^{(1)}&\ldots&a_{1}^{(k)}\\\
\vdots&&\\\
a_{n}^{(1)}&\ldots&a_{n}^{(k)}\end{array}\right)\biggl{|}\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right).\end{split}$

Here $\eta_{j}(\theta_{l})$ are the independent Poisson random variables with
mean $\theta_{l}/j$, i.e.

(6.4)
$\operatorname{Prob}\left\\{\eta_{j}(\theta_{l})=p\right\\}=e^{-\theta_{l}/j}\left(\frac{\theta_{l}}{j}\right)^{p}\frac{1}{p!},\;\;p=0,1,2,\ldots.$

###### Proof.

The fact that (1.3) defines a probability distribution on the set of all
matrices (6.1) satisfying (6.2) implies

(6.5)
$\sum\limits_{b_{j}^{(l)}:\;\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}jb_{j}^{(l)}=n}\prod_{l=1}^{k}\prod\limits_{j=1}^{n}\frac{\left(\theta_{l}/j\right)^{b_{j}^{(l)}}}{(b_{j}^{(l)})!}=\frac{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}{n!}.$

Note that the left-hand side of (6.5) can be understood as

$\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}.$

We will show that

(6.6)
$\begin{split}&\operatorname{Prob}\left\\{A_{j}^{(l)}(n)=a_{j}^{(l)};\;j=1,\ldots,n;\;l=1,\ldots,k\right\\}\\\
&=\operatorname{Prob}\left\\{\eta_{j}(\theta_{l})=a_{j}^{(l)};\;j=1,\ldots,n;\;l=1,\ldots,k\biggl{|}\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}.\end{split}$

Indeed, the right-hand side of equation (6.6) can be computed as

(6.7)
$\begin{split}&\operatorname{Prob}\left\\{\eta_{j}(\theta_{l})=a_{j}^{(l)};\;j=1,\ldots,n;\;l=1,\ldots,k\biggl{|}\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}\\\
&=\frac{\operatorname{Prob}\left\\{\eta_{j}(\theta_{l})=a_{j}^{(l)};\;j=1,\ldots,n;\;l=1,\ldots,k\right\\}}{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}}\\\
&=\frac{n!}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\prod\limits_{l=1}^{k}\prod\limits_{j=1}^{n}\frac{\left(\theta_{l}/j\right)^{a_{j}^{(l)}}}{\left(a_{j}^{(l)}\right)!}e^{-\theta_{l}/j},\end{split}$

which is equal to the probability in the left-hand side of equation (6.6). ∎

Our next result is Theorem 6.2 which says that as $n\rightarrow\infty$ the
random matrix $A(n)$ can be approximated by a matrix whose entries are
independent Poisson random variables.

######  Theorem 6.2.

For any fixed strictly positive integer $m$, it holds that

(6.8) $\left(\begin{array}[]{ccc}A_{1}^{(1)}(n)&\ldots&A_{1}^{(k)}(n)\\\
\vdots&&\\\
A_{m}^{(1)}(n)&\ldots&A_{m}^{(k)}(n)\end{array}\right)\overset{D}{\longrightarrow}\left(\begin{array}[]{ccc}\eta_{1}(\theta_{1})&\ldots&\eta_{1}(\theta_{k})\\\
\vdots&&\\\
\eta_{m}(\theta_{1})&\ldots&\eta_{m}(\theta_{k})\end{array}\right),$

as $n\rightarrow\infty$. Here $\eta_{j}(\theta_{l})$ are the independent
Poisson random variables with mean $\theta_{l}/j$.

###### Proof.

We want to compute the asymptotics of

(6.9)
$\operatorname{Prob}\left\\{A_{j}^{(l)}(n,\theta_{l})=a_{j}^{(l)},\;j=1,\ldots,m,l=1,\ldots,k\right\\}$

as $n\rightarrow\infty$, where $a_{j}^{(l)}$ are non-negative integers
satisfying condition (6.2). By Proposition 6.1 we can rewrite (6.9) as

(6.10)
$\begin{split}&\operatorname{Prob}\left\\{\eta_{j}(\theta_{l})=a_{j}^{(l)},\;j=1,\ldots,m,\;l=1,\ldots,k\right\\}\\\
&\times\frac{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}j\eta_{j}(\theta_{l})=n-a^{(1)}-\ldots-a^{(k)}\right\\}}{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}},\end{split}$

where $a^{(l)}=\sum\limits_{j=1}^{m}ja_{j}^{(l)}$, $l=1,\ldots,k$. It is
enough to show that

(6.11)
$\begin{split}&\frac{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}j\eta_{j}(\theta_{l})=n-a^{(1)}-\ldots-a^{(k)}\right\\}}{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}}\longrightarrow
1,\end{split}$

as $n\rightarrow\infty$. Denote

(6.12)
$T_{0,n}(\theta_{l})=\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})\;\;T_{m,n}(\theta_{l})=\sum\limits_{j=m+1}^{n}j\eta_{j}(\theta_{l}).$

The right-hand side of (6.10) has the form

(6.13)
$\frac{\operatorname{Prob}\left\\{T_{m,n}(\theta_{1})+\ldots+T_{m,n}(\theta_{k})=n-a^{(1)}-\ldots-a^{(k)}\right\\}}{\operatorname{Prob}\left\\{T_{0,n}(\theta_{1})+\ldots+T_{0,n}(\theta_{k})=n\right\\}}.$

Formula (6.5) implies

(6.14)
$\operatorname{Prob}\left\\{T_{0,n}(\theta_{1})+\ldots+T_{0,n}(\theta_{k})=n\right\\}=\frac{(\theta_{1}+\ldots+\theta_{k})_{n}}{n!}e^{-(\theta_{1}+\ldots+\theta_{k})\sum_{j=1}^{n}j^{-1}}.$

In order to proceed we wish to compute the generating function of the random
variable $\sum\limits_{l=1}^{k}T_{m,n}(\theta_{l})$. Using equation (6.12) we
obtain

(6.15)
$\begin{split}&\mathbb{E}\left(x^{\sum\limits_{l=1}^{k}T_{m,n}(\theta_{l})}\right)=\mathbb{E}\left(x^{\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}j\eta_{j}(\theta_{l})}\right)=\prod\limits_{l=1}^{k}\prod\limits_{j=m+1}^{n}\mathbb{E}\left(x^{j\eta_{j}(\theta_{l})}\right)\\\
&=\prod\limits_{l=1}^{k}\prod\limits_{j=m+1}^{n}e^{\frac{\theta_{l}}{j}(x^{j}-1)}.\end{split}$

For a positive integer-valued random variable $Z$ the following formula holds
true

(6.16)
$\operatorname{Prob}\left(Z=k\right)=\frac{1}{k!}\left(\frac{d^{k}}{dx^{k}}\mathbb{E}\left(x^{Z}\right)\right)\biggl{|}_{x=0}.$

Using (6.15) and (6.16) we obtain the formula

(6.17)
$\begin{split}&(n-a^{(1)}-\ldots-a^{(k)})!\operatorname{Prob}\left(\sum\limits_{l=1}^{k}T_{m,n}(\theta_{l})=n-a^{(1)}-\ldots-a^{(k)}\right)\\\
&=e^{-\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}\frac{\theta_{l}}{j}}\biggl{[}\frac{d^{n-a^{(1)}-\ldots-a^{(k)}}}{dx^{n-a^{(1)}-\ldots-a^{(k)}}}e^{\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}\frac{\theta_{l}x^{j}}{j}}\biggr{]}\biggr{|}_{x=0}.\end{split}$

Set

(6.18)
$\begin{split}&g_{l}(x)=e^{-\theta_{l}\sum\limits_{j=1}^{m}\frac{x^{j}}{j}},\;l=1,\ldots,k;\;\;G(x)=g_{1}(x)\ldots
g_{k}(x).\end{split}$

Then

(6.19)
$\begin{split}&\operatorname{Prob}\left(\sum\limits_{l=1}^{k}T_{m,n}(\theta_{l})=n-a^{(1)}-\ldots-a^{(k)}\right)\\\
&=\frac{e^{-\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}\frac{\theta_{l}}{j}}}{(n-a^{(1)}-\ldots-a^{(k)})!}\biggl{[}\frac{d^{n-a^{(1)}-\ldots-a^{(k)}}}{dx^{n-a^{(1)}-\ldots-a^{(k)}}}\left(G(x)e^{-\left(\theta_{1}+\ldots+\theta_{k}\right)\log(1-x)}\right)\biggr{]}\biggr{|}_{x=0}.\end{split}$

If $A=a^{(1)}+\ldots+a^{(k)}$, and $\Omega=\theta_{1}+\ldots+\theta_{k}$, then

$\begin{split}&\biggl{[}\frac{d^{n-A}}{dx^{n-A}}\left(G(x)e^{-\Omega\log(1-x)}\right)\biggr{]}\biggl{|}_{x=0}=G(1)\left(\Omega\right)_{n-A}-G^{\prime}(1)(\Omega-1)_{n-A}+\frac{G^{{}^{\prime\prime}}(1)}{2}(\Omega-2)_{n-A}-\ldots\\\
&=e^{-\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{m}\theta_{l}}\left(\Omega\right)_{n-A}\left[1+\frac{\left(m\sum\limits_{l=1}^{k}\theta_{l}\right)(\Omega-1)}{\Omega-1+n-A}+O\left(\frac{1}{n^{2}}\right)\right].\end{split}$

Therefore,

(6.20)
$\begin{split}&\operatorname{Prob}\left(\sum\limits_{l=1}^{k}T_{m,n}(\theta_{l})=n-a^{(1)}-\ldots-a^{(k)}\right)=\frac{e^{-\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}\frac{\theta_{l}}{j}}\left(\theta_{1}+\ldots+\theta_{k}\right)_{n-a^{(1)}-\ldots-a^{(k)}}}{(n-a^{(1)}-\ldots-a^{(k)})!}\\\
&\times\left[1+\frac{\left(m\sum\limits_{l=1}^{k}\theta_{l}\right)(\theta_{1}+\ldots+\theta_{k}-1)}{\theta_{1}+\ldots+\theta_{k}-1+n-a^{(1)}-\ldots-a^{(k)}}+O\left(\frac{1}{n^{2}}\right)\right].\end{split}$

We conclude that

(6.21)
$\begin{split}&\frac{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=m+1}^{n}j\eta_{j}(\theta_{l})=n-a^{(1)}-\ldots-a^{(k)}\right\\}}{\operatorname{Prob}\left\\{\sum\limits_{l=1}^{k}\sum\limits_{j=1}^{n}j\eta_{j}(\theta_{l})=n\right\\}}\\\
&=\frac{n!}{(n-a^{(1)}-\ldots-a^{(k)})!}\;\frac{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n-a^{(1)}-\ldots-a^{(k)}}}{\left(\theta_{1}+\ldots+\theta_{k}\right)_{n}}\\\
&\times\left[1+\frac{\left(m\sum\limits_{l=1}^{k}\theta_{l}\right)(\theta_{1}+\ldots+\theta_{k}-1)}{\theta_{1}+\ldots+\theta_{k}-1+n-a^{(1)}-\ldots-a^{(k)}}+O\left(\frac{1}{n^{2}}\right)\right]\longrightarrow
1,\end{split}$

as $n\rightarrow\infty$. Theorem 6.2 is proved. ∎

######  Remark 6.3.

In the situation where the mutation parameter grows with the sample size $n$
the Poisson approximation for the Ewens sampling formula was considered in
Tsukuda [22]. Methods of Ref. [22] can be applied to the refined Ewens
sampling formula (equation (1.3)) as well to study asymptotic regime where the
mutation parameters $\theta_{1}$, $\ldots$, $\theta_{k}$ grow with $n$ from
the viewpoint of Poisson approximation. We postpone this issue to a future
work.

## References

  * [1] Aldous, D. J. Exchangeability and related topics. $\acute{\mbox{E}}$cole d’$\acute{\mbox{e}}$t$\acute{\mbox{e}}$ de probabilit$\acute{\mbox{e}}$s de Saint-Flour, XIII—1983, 1–198, Lecture Notes in Math., 1117, Springer, Berlin, 1985.
  * [2] Arratia, R.; Barbour, A.; Tavar$\acute{\mbox{e}}$, S. Poisson process approximation for the ewns sampling formula. Ann. Appl. Probab. 2 (1992), no. 3, 519-–535.
  * [3] Arratia, R.; Barbour, A.; Tavar$\acute{\mbox{e}}$, S. Logarithmic combinatorial structures: a probabilistic approach. European Mathematical Society, Zurich, 2003.
  * [4] Crane, H. The Ubiquitous Ewens Sampling Formula. Statistical Science 31 (2016) 1–19.
  * [5] Durrett, R. Probability Models for DNA Sequence Evolution. Probability and Its Applications. Springer 2008.
  * [6] Etheridge, A. Some Mathematical Models from Population Genetics. Lecture Notes in Mathematics. Springer 2012.
  * [7] Ewens, W. J. The sampling theory of selectively neutral alleles. Theoret. Population Biol. 3 (1972).
  * [8] Griffiths, R. C.; Lesard, S. Ewens’ sampling formula and related formulae: combinatorial proofs, extension to a variable population size and applications to ages of alleles. Theor. Population Biology 68 (2005) 167–177.
  * [9] Hoppe, F. M. Pólya-like urns and the Ewens’ sampling formula. J. Math. Biol. 20 (1984), no. 1, 91–-94.
  * [10] Hoppe, F. M. The sampling theory of neutral alleles and an urn model in population genetics. J. Math. Biol. 25 (1987), no. 2, 123–-159.
  * [11] Feng, S. Large deviations associated with Poisson-Dirichlet distribution and Ewens sampling formula. Ann. Appl. Probab. 17 (2007), no. 5–6, 1570–-1595.
  * [12] Feng, S. The Poisson-Dirichlet Distribution and Related Topics. Models and Asymptotic Behaviors. Springer 2010.
  * [13] Kingman, J. F. C. Random partitions in population genetics. Proc. R. Soc. London A 361, (1978), 1–20.
  * [14] Kingman, J. F. C. The representation of partition structures. J. London Math. Soc. 18 (1978), 374–380.
  * [15] Macdonald, I. Symmetric Functions and Hall Polynomials. Second edition, Oxford Mathematical Monographs, Oxford University Press (1995).
  * [16] Mahmoud, H. M. P$\acute{\mbox{o}}$lya Urn Models. Texts in Statistical Sciences (2009).
  * [17] O’Reilly, R.; Elphick, H. E. Development, clinical utility, and place of ivacaftor in the treatment of cystic fibrosis. Drag Design, Development and Therapy 7 (2013), 929–937.
  * [18] Strahov, E. Multiple partition structures and harmonic functions on branching graphs. Adv. in Appl. Math. 153 (2024), Paper No. 102617, 49 pp.
  * [19] Strahov, E. Generalized regular representations of big wreath products. Israel J. of Math. (2024) (will appear)
  * [20] Tavar$\acute{\mbox{e}}$, S. The magical Ewens sampling formula. Bull. London Math. Soc. 53 (2021), 1563–1582.
  * [21] Tsukuda, K. Estimating the large mutation parameter of the Ewens sampling formula. J. Appl. Probability 54 (2017), 42–54.
  * [22] Tsukuda, K. On Poisson approximations for the Ewens sampling formula when the mutation parameter grows with the sample size. Ann. Appl. Probab. 29 (2019), no. 2, 1188–1232.
  * [23] Watterson, G. A. The stationary distribution of the infinitely-many neutral alleles diffusion model. J. Appl. Probability 13 (1976), no. 4, 639–-651.