# An evolutionary model that satisfies detailed balance

Jüri Lember111Institute of Mathematics and Statistics, University of Tartu,
Estonia, Chris Watkins222Department of Computer Science, Royal Holloway,
University of London, UK

###### Abstract

We propose a class of evolutionary models that involves an arbitrary
exchangeable process as the breeding process and different selection schemes.
In those models, a new genome is born according to the breeding process, and
then a genome is removed according to the selection scheme that involves
fitness. Thus the population size remains constant. The process evolves
according to a Markov chain, and, unlike in many other existing models, the
stationary distribution – so called mutation-selection equilibrium – can be
easily found and studied. The behaviour of the stationary distribution when
the population size increases is our main object of interest. Several phase-
transition theorems are proved.

##### Keywords:

Markov chain Monte Carlo, Dirichlet distribution, de Finetti theorem, weak
convergence of probability measures.

##### AMS classification:

60J10, 60B10

##  1 Introduction

We introduce a probability model of evolution that has three purposes: as an
abstract model of biological evolution; as a class of efficiently
implementable genetic algorithms that are easy to analyse theoretically; and
as a link between genetic algorithms and Bayesian non-parametric MCMC methods.
The stationary distribution of the model can be expressed in closed form for
arbitrary fitness functions: this enables us to investigate the behaviour of
the model for different population sizes, mutation rates, and fitness
scalings. We find two phase transitions that occur for all fitness functions.
Our approach is most applicable to evolution by sexual rather than asexual
reproduction.

In any model of evolution under constant conditions, there is a temporal
sequence of possibly overlapping populations. In the transition from each
population to the next, one or more new individuals are ‘born’, and the same
number of individuals are removed from the population, or ‘die’. Each new
generation is a population that depends only upon the previous parent
population, so that the sequence of populations is a Markov chain. In a model
with mutation, the Markov chain is irreducible and has a unique stationary
distribution that is also known as the mutation-selection equilibrium.
Typically we wish to know what the mutation-selection equilibrium is.
Unfortunately the mutation-selection equilibrium is notoriously hard to
characterize, even for apparently simple and elegant models of breeding,
mutation, and selection. The reason is that in previous models of evolution
with sexual reproduction, the Markov chain of populations is irreversible, so
that there is no obvious method of finding the stationary distribution other
than attempting to compute eigenvectors of the transition matrix directly, as
done in [18], but these calculations are neither easy nor revealing for
arbitrary fitness functions.

For a reversible Markov chain, the stationary distribution may be found by
verifying detailed balance conditions. In our model, introduced in [19], we
start by writing the mutation-selection equilibrium in a convenient closed
form, and then exhibit MCMC kernels that implement reversible Markov chains
with this stationary distribution, and for which the proposal and acceptance
algorithms are plausible abstract models of breeding, mutation, and selection.

Each generation starts with a population of $n$ individuals. One new
individual is ‘bred’ – that is, sampled conditionally on the existing
population – to produce an expanded population of $n+1$ individuals; from this
expanded population, one individual is selected to be discarded, leaving a new
population of size $n$ to start the next generation. This might appear similar
to the Moran Process [12] but in fact it is very different, as explained in
section 3.6

In previous evolutionary models such as [4], and in genetic algorithms such as
[11], the reproductive fitness of a genome is modelled as the genome’s rate of
breeding, and not as its probability of death. In the breeding phase of each
generation, fit genomes are chosen to breed more frequently than unfit
genomes: in these models, discarding individuals – or ‘death’ – is modelled as
random deletion from the population.

Here we model breeding as conditional sampling in which all existing members
of the population are treated equally, regardless of their fitness. We discard
individuals according to their fitness: less-fit individuals are more likely
to be discarded, so that they remain in the population for a shorter time, and
they are therefore sampled less as a result of their shorter lifetime, and so
contribute less to the new individuals that are ‘bred’. Thus differences in
reproductive fitness are modelled as differences in longevity. This is a
significant design choice in our models because breeding and selection are
modelled separately, and this turns out to greatly simplify the analysis.

As explained in detail in sections 2 and 3, careful choices in modelling
breeding as conditional sampling, and in modelling fitness by stochastic
rejection of the less fit, enable us to construct a reversible Markov chain of
populations that is a form of Metropolis-Hastings process.

Throughout the paper, we suppose there is a finite set ${\cal X}$ of possible
genomes, and $\xi_{1},\xi_{2},\ldots$ denotes an exchangeable ${\cal
X}$-valued stochastic process. For any population of $n$ genomes
$\mathbf{x}=(x_{1},\ldots,x_{n})\in{\cal X}^{n}$, we define

$P(\xi_{1}=x_{1},\ldots,\xi_{n}=x_{n})=:P_{\xi}(x_{1},\ldots,x_{n})=:P_{\xi}(\mathbf{x})$
(1.1)

By definition of exchangeability, we have that for any permutation $\sigma$,

$P_{\xi}(x_{1},\ldots,x_{n})=P_{\xi}(x_{\sigma(1)},\ldots,x_{\sigma(n)}).$

Given a population of genomes $\mathbf{x}=(x_{1},\ldots,x_{n})$, we breed an
$n+1$’th genome by sampling $x_{n+1}$ conditionally:

$x_{n+1}\sim P_{\xi}(\cdot\mid x_{1},\ldots,x_{n}).$

The ‘fitness’ of a genome $x$ is denoted by $w(x)>0$, where $w$ is an
arbitrary strictly positive function over $\mathcal{X}$. In the evolutionary
models we propose, in each generation one individual is ‘bred’ by conditional
sampling from the existing population, and then an individual is discarded in
a fitness-biased way, so that less-fit individuals are more likely to be
discarded. The stationary distribution of populations factorises into the
form:

$\underbrace{P_{n}(x_{1},\ldots,x_{n})}_{\text{stationary
distribution}}=\frac{1}{Z_{n}}\underbrace{P_{\xi}(x_{1},\ldots,x_{n})}_{\text{breeding
term}}~{}\underbrace{w(x_{1})\cdots w(x_{n})}_{\text{fitness term}}.$ (1.2)

The process is similar to, but not the same as non-parametric Bayesian MCMC,
and we make this connection explicit in section 3.5.

A genetic algorithm has three basic parameters: population size, mutation
rate, and fitness scaling. The most important results of this paper
characterise the interacting effect of these parameters as population size
$n\rightarrow\infty$.

We particularly consider exchangeable processes $\xi$ that are (collections
of) Polya urn processes, also known as Dirichlet-categorical processes because
these admit a notion of ‘mutation’. Conditional sampling from these processes
is directly interpretable as a simplified model of sexual breeding with
mutation, as explained in section 2. We denote the parameters defining the
Dirichlet prior(s) as $\alpha$. Since fitness $w>0$, we often write $w$ in
terms of a function $\phi:\mathcal{X}\rightarrow\mathbb{R}$ such that
$w(x)=\exp(-\phi(x))$. We study the interaction of $n$, $\alpha$, and $\phi$
on the stationary distribution by letting population size $n\rightarrow\infty$
while the fitness of each genome $x$ scales with $n$ as
$\exp(n^{-\lambda}\phi(x))$, and $\alpha$ scales with $n$ as
$n^{1-\lambda}\alpha$; in sections 5 and 6 we use the de Finetti
representation of $\xi$ to derive limiting forms of the marginal distributions
over $\mathcal{X}$ of

$P_{n}(x_{1},\ldots,x_{n}\mid
n^{1-\lambda}\alpha;n^{-\lambda}\phi)~{}~{}\text{ as $n\rightarrow\infty$}.$

We establish that with this rescaling, for any fitness function $\phi$ over
$\mathcal{X}$, and for any Dirichlet-categorical process $\xi$, there are
three different non-trivial population distributions in the limit as
$n\rightarrow\infty$:

Constant mutation-rate limit

: $\lambda=0$, so that $\alpha$ scales as $n\alpha$ and $\phi$ is unscaled;

Low mutation, weak selection limit

: $\lambda\in(0,1]$, so that $\alpha$ scales as $n^{1-\lambda}\alpha$ and
$\phi$ as $n^{-\lambda}\phi$;

Critical low mutation limit

: $\lambda=1$, so that $\alpha$ is unscaled, and $\phi$ scales as
$\frac{\phi}{n}$.

The phase transitions at $\lambda=0$ and $\lambda=1$ are sharp. We
characterise the limiting distributions in the case that there is a unique
maximal element of the posterior distribution. As far as we are aware, these
are the first results that give exact explicit expressions for stationary
distributions of genetic algorithms with arbitrary fitness functions.

The paper is organised as follows. In section 2 we give examples of
exchangeable conditional sampling procedures that can be regarded as abstract
models of breeding, and section 3 gives examples of selection procedures that,
together with any of the breeding procedures, will produce a reversible Markov
chain of populations with the stationary distribution of equation 1.2. Section
3.4 establishes that the stationary distributions are invariant to
multiplicative fitness noise; section 3.5 shows that special cases of this
process are forms of Bayesian inference by MCMC. This family of evolutionary
algorithms might appear similar to the Moran Process – but section 3.6
explains that the Moran process is quite different. Next, section 4
establishes basic conditions on the convergence of the stationary distribution
to a limit distribution as population size tends to infinity. Sections 5 and 6
characterise the limiting forms of the stationary distribution for large
populations. These are the main results of this paper. We present
computational experiments demonstrating our results in section 7. Finally we
discuss implications of our results for genetic algorithms and evolutionary
modelling in 8.

##  2 Breeding and mutation

We model breeding as Gibbs sampling [7] from an exchangeable distribution;
exchangeable Gibbs sampling is a standard technique in statistical
nonparametric MCMC methods, for example [14, 10], but in that context it is
not of course regarded as a model of breeding. The property of exchangeability
of $P_{\xi}$ will be used in two ways: first, in section 3 we will use it to
establish detailed balance for several selection procedures, which establishes
that the stationary distribution of populations is indeed as given in equation
1.2. Second, in section 4 we use the de Finetti integral representation of
$\xi$ to establish limit properties of the stationary distribution as
$n\rightarrow\infty$.

We now give examples of conditional sampling that can be regarded as plausible
models of breeding.

##### Dirichlet-Categorical Process.

In the simplest case, each ‘genome’ consists of only one ‘gene’ which can be
one of $K$ possible alleles, that we denote by $\\{1,\ldots,K\\}$. The
exchangeable process $\xi$ is the well known Polya urn model for a Dirichlet
process with discrete base distribution. We recall the definition of this
process. Let $\alpha=(\alpha_{1},\ldots,\alpha_{K})$, $\alpha_{i}>0$ be the
prior parameters of the base distribution; we write
$|\alpha|:=\alpha_{1}+\cdots+\alpha_{K}$. Let $\xi_{1},\xi_{2},\ldots$ be a
random process over $\\{1,\ldots,K\\}$. Let $n\geq 0$. Given
$\mathbf{x}=x_{1},\ldots,x_{n}$, we denote the number of $k$-s in the sequence
by $n_{k}(\mathbf{x})$:

$n_{k}=\sum_{i=1}^{n}I_{\\{k\\}}(x_{i}).$ (2.1)

Thus

$P_{\xi}(\xi_{n+1}=k\mid
x_{1},\ldots,x_{n},\alpha):=\frac{n_{k}+\alpha_{k}}{n+|\alpha|}.$

It follows that

$P_{\xi}(x_{1},\ldots,x_{n}\mid\mathbf{\alpha})=\frac{\alpha_{1}(\alpha_{1}+1)\cdots(\alpha_{1}+n_{1}-1)\cdots\alpha_{K}(\alpha_{K}+1)\cdots(\alpha_{K}+n_{K}-1)}{|\alpha|(|\alpha|+1)\cdots(|\alpha|+n-1)}$

and $\xi_{1},\xi_{2},\ldots$ is infinitely exchangeable by inspection. By de
Finetti’s theorem

$P_{\xi}(x_{1},\ldots,x_{n}\mid\mathbf{\alpha})=\int_{{\cal P}}q(x_{1})\cdots
q(x_{n})\pi(dq),$

where ${\cal P}=\\{(q_{1},\ldots,q_{K}):q_{i}\geq 0,\sum q_{i}=1\\}$ is the
set of all probability vectors (simplex) and the prior measure $\pi$ is
Dirichlet distribution, i.e. $\pi={\rm Dir}(\alpha_{1},\ldots,\alpha_{K})$.
This process is in a sense the central object of our study.

##### Concentration parameter and mutation rate.

The concentration parameter ${|\alpha|}=\alpha_{1}+\cdots+\alpha_{K}$ may be
viewed as determining a mutation rate that depends on $n$. With $n$ balls in
the urn, when a new ball is sampled, there is a probability
$\frac{|\alpha|}{n+|\alpha|}$ that the ball will be sampled from the
collection of ‘prior’ balls, rather than the actual balls in the urn. This
probability is independent of the colors of $n$ actual balls present: since a
new colour may be introduced in this way, we regard this as analogous to a
mutation. The mutation rate $u$ is

$u=\frac{|\alpha|}{|\alpha|+n}~{}~{}~{}\text{or
equivalently}~{}~{}~{}|\alpha|=\frac{nu}{1-u}$

In our evolutionary processes, one new genome is sampled, and one genome is
then discarded at each generation, so that $n$ remains constant. To define
processes with the same mutation rate and different values of $n$, $\alpha$
must be adjusted to depend on $n$.

Note that in this model, mutations only occur at birth, and – importantly –
the distribution of mutations does not depend on the frequencies of different
colours that are currently in the population; mutations are distributed
according to a fixed prior distribution determined by the frequencies of
‘magic’ balls in the urn. This is strictly less general than an alternative
model in which breeding occurs as follows: a ball $x$ is sampled from the urn,
and then a mutated ball $x^{\dagger}\sim P_{M}(\cdot\mid x)$ is conditionally
sampled according to some mutation distribution $P_{M}$, and then $x$ and the
new mutated ball $x^{\dagger}$ are both returned to the urn. But this
alternative model would not necessarily be reversible, and we do not consider
it further.

In our model, we count as a mutation any ball which results from a draw of a
‘magic’ ball: this definition is consistent with an extension of our model to
Dirichlet processes with continuous base distributions, which we intend to
consider in future work.

###  2.1 Complex genomes: direct product of Dirichlet processes

More complex evolutionary models such as genetic algorithms require more
complex genomes. Suppose that each genome is a vector of $L$ genes,
$x_{i}=(y^{1}_{i},\ldots,y^{L}_{i})$. Let $\xi$ be a direct product of $L$
independent exchangeable processes $\xi=(\xi^{1},\ldots,\xi^{L})$. Then

$P_{\xi}(x_{1},\ldots,x_{n}):=\prod_{j=1}^{L}P_{\xi^{j}}(y^{j}_{1},\ldots,y^{j}_{n})$
(2.2)

which clearly is exchangeable as well. Exchangeable sampling from
$P_{\xi}(\cdot\mid x_{1},\ldots,x_{n})$, where each $x_{i}$ is a vector of
discrete values, can be viewed as a model of sexual reproduction with the
assumption of linkage equilibrium. A new vector
$(y_{n+1}^{1},\ldots,y_{n+1}^{L})$ is sampled by, for $1\leq j\leq L$,
sampling the $j$-th component $y_{n+1}^{j}\sim P_{\xi^{j}}(\cdot\mid
y^{j}_{1},\ldots,y^{j}_{n})$ independently from the rest of the components. In
words, each new element of the vector $x_{n+1}$ is either a copy of the
corresponding element of a randomly chosen member of the existing population,
or else a mutation. Instead of a new ‘child’ genome being constructed by
random recombination of two parent genomes, it is instead a random
recombination of all $n$ existing genomes in the population, with mutations.

This method of constructing new genomes by ‘$n$-way recombination’ is a widely
used approach in genetic algorithms, as used by [1, 2, 3] and others, and
sexual reproduction with full linkage equilibrium is a standard simplified
model of sexual reproduction in population genetics theory [4, 6].

The extension to Cartesian products of Dirichlet processes might appear rather
simple because each component of a new genome is sampled independently of the
others; however, this extension can lead to models of great complexity because
the fitness function $w$, or equivalently $\phi$, can be an arbitrary function
on $\mathcal{X}^{L}$, so that the stationary distribution need not be a
product distribution. In genetic language, the fitness function can have
arbitrary epistasis.

Note that there are other discrete exchangeable distributions based on
Dirichlet distributions that could also be used as $P_{\xi}$. A notable
example is the discrete fragmentation-coagulation sequence process introduced
in [5]; this was intended as statistical model for imputing phasing in genetic
analysis, but it could also be used as a breeding distribution for our
purposes.

##  3 Selection

Several MCMC sampling methods give the factorized stationary distribution
given by equation (1.2) and at the same time are models of sexual reproduction
that are as plausible as those used in evolutionary computation or simplified
models in population genetics.

We suppose that each element $x\in{\cal X}$ has a strictly positive weight
$w(x)$. In context, we will denote the weights $w(x_{1}),\ldots,w(x_{n})$ as
$w_{1},\ldots,w_{n}$.

###  3.1 Single tournament selection

We suppose that when a new genome $x_{n+1}$ is ‘born’ and added to the
population, it competes to survive by having a tournament with another
randomly selected member of the population, $x_{i}$ say. The probability that
$x_{n+1}$ wins the tournament and ejects $x_{i}$ from the population is
$\frac{w(x_{n+1})}{w(x_{n+1})+w(x_{i})}$. This probability is always well
defined since $w$ is strictly positive. An equally valid tournament winning
probability is that $x_{n+1}$ wins with probability
$\min\\{1,\frac{w(x_{n+1})}{w(x_{i})}\\}$. These two tournament winning
probabilities are simply different formulations of the Metropolis-Hastings
acceptance rule. The proof below establishes detailed balance for the first
winning rule. The algorithm for performing one generation of breeding,
mutation, and selection is:

  1. 1.

Sample $x_{n+1}\sim P_{\xi}(\cdot\mid x_{1},\ldots,x_{n})$

  2. 2.

Sample $i$ randomly from $\\{1,\ldots,n\\}$

  3. 3.

With probability $\frac{w_{n+1}}{w_{i}+w_{n+1}}$ replace $x_{i}$ with
$x_{n+1}$ and discard $x_{i}$, otherwise discard $x_{n+1}$.

Let $\mathbf{x}$, $\mathbf{x}^{\prime}$ be populations defined as

$\displaystyle\mathbf{x}$
$\displaystyle=x_{1},\ldots,x_{n}\quad\mathbf{x}^{\prime}=x_{1},\ldots,x_{i-1},x_{n+1},x_{i+1},\ldots,x_{n}$

Recall the measure $P_{n}(\mathbf{x})$ defined in (1.2). We now show that this
measure satisfies detailed balance. By exchangeability of $\xi$, we have:

$P_{\xi}(\mathbf{x})P_{\xi}(x_{n+1}\mid\mathbf{x})=P_{\xi}(x_{1},\ldots,x_{n+1})=P_{\xi}(\mathbf{x}^{\prime})P_{\xi}(x_{i}\mid\mathbf{x}^{\prime})$
(3.1)

Note that $x_{n+1}$ has tournament against $x_{i}$ with probability
$\frac{1}{n}$ and wins with probability $\frac{w_{n+1}}{w_{i}+w_{n+1}}$, so:

$\displaystyle P_{n}(\mathbf{x})P(\mathbf{x}\to\mathbf{x}^{\prime})$
$\displaystyle=P_{n}(\mathbf{x})\cdot
P_{\xi}(x_{n+1}\mid\mathbf{x})\cdot\frac{1}{n}\frac{w_{n+1}}{w_{i}+w_{n+1}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(\mathbf{x})w_{1}\cdots w_{n}\cdot
P_{\xi}(x_{n+1}\mid\mathbf{x})\cdot\frac{1}{n}\frac{w_{n+1}}{w_{i}+w_{n+1}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(x_{1},\ldots,x_{n+1})\cdot\frac{1}{n}\frac{w_{1}\cdots
w_{n+1}}{w_{i}+w_{n+1}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(\mathbf{x}^{\prime})P_{\xi}(x_{i}\mid\mathbf{x}^{\prime})\cdot\frac{1}{n}\frac{w_{1}\cdots
w_{n+1}}{w_{i}+w_{n+1}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(\mathbf{x}^{\prime})w_{1}\cdots
w_{i-1}w_{n+1}w_{i+1}\cdots w_{n}\cdot
P_{\xi}(x_{i}\mid\mathbf{x}^{\prime})\cdot\frac{1}{n}\frac{w_{i}}{w_{i}+w_{n+1}}$
$\displaystyle=P_{n}(\mathbf{x}^{\prime})P(\mathbf{x}^{\prime}\to\mathbf{x}).$

###  3.2 Inverse fitness selection: limit of many tournaments

Suppose that many tournaments are fought, and each time the loser of the
previous tournament fights another randomly chosen genome from the population.
After many tournaments, and at a stopping time, the current loser is ejected.
The current loser evolves according to a irreducible aperiodic Markov chain,
and the limiting distribution of ejection is the stationary distribution of
that chain:

$P(\text{ eject $i$
})=\frac{\frac{1}{w_{i}}}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}.$ (3.2)

The algorithm for performing one generation of breeding and selection is then:

  1. 1.

Sample $x_{n+1}\sim P_{\xi}(\cdot\mid x_{1},\ldots,x_{n})$

  2. 2.

Sample $i$ from a discrete p.d. over $\\{1,\ldots,n+1\\}$ with probabilities
proportional to $\\{\frac{1}{w_{1}},\ldots,\frac{1}{w_{n+1}}\\}$

  3. 3.

Discard $x_{i}$

This process too satisfies detailed balance. With $\mathbf{x}$ and
$\mathbf{x}^{\prime}$ defined as above, note that:

$\displaystyle P_{n}(\mathbf{x})P(\mathbf{x}\to\mathbf{x}^{\prime})$
$\displaystyle=P_{n}(\mathbf{x})\cdot
P_{\xi}(x_{n+1}\mid\mathbf{x})\cdot\frac{\frac{1}{w_{i}}}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(\mathbf{x})w_{1}\cdots w_{n}\cdot
P_{\xi}(x_{n+1}\mid\mathbf{x})\cdot\frac{\frac{1}{w_{i}}}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(x_{1},\ldots,x_{n+1})\frac{w_{1}\cdots
w_{n+1}}{w_{n+1}}\frac{\frac{1}{w_{i}}}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}$
$\displaystyle=\frac{1}{Z_{n}}P_{\xi}(x_{1},\ldots,x_{n+1})\frac{w_{1}\cdots
w_{n+1}}{w_{i}}\frac{\frac{1}{w_{n+1}}}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}$
$\displaystyle=P_{n}(\mathbf{x}^{\prime})\cdot
P_{\xi}(x_{i}\mid\mathbf{x}^{\prime})\cdot\frac{\frac{1}{w_{n+1}}}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}$
$\displaystyle=P_{n}(\mathbf{x}^{\prime})P(\mathbf{x}^{\prime}\to\mathbf{x}).$

Note in passing that

$\mathbf{E}\\{\mbox{weight of rejected genome}\mid
w_{1},\ldots,w_{n+1}\\}=\frac{1}{\frac{1}{w_{1}}+\cdots+\frac{1}{w_{n+1}}}.$
(3.3)

That is, the expected weight of the rejected genome is the harmonic mean of
the weights of the genomes in the current population, including the newly
added genome $x_{n+1}$.

##### The case $n=1$ reduces to Metropolis-Hastings:

with a population of size 1, we have proposal distribution
$P_{\xi}(x^{\prime}\mid x)P_{\xi}(x)=P_{\xi}(x,x^{\prime})=P_{\xi}(x\mid
x^{\prime})P_{\xi}(x^{\prime})$ and acceptance probability of $x^{\prime}$
given $x$ of $\frac{w(x^{\prime})}{w(x)+w(x^{\prime})}$, which gives a
stationary distribution

$P_{1}(x)=\frac{P_{\xi}(x)w(x)}{\sum_{x}P_{\xi}(x)w(x)}.$

###  3.3 Exchangeable breeding of many offspring

Another MCMC process which is also interpretable as an evolutionary algorithm
breeds some arbitrary number $m$ of offspring to give a population of size
$n+m$, and then from these selects $n$ genomes to form the next generation.
The algorithm for a single generation is as follows:

  1. 1.

Pick a random number $m$ of offspring to breed, and a random number $t$ of
tournaments to conduct. Both $m$ and $t$ should be independent of the current
population $x_{1},\ldots,x_{n}$

  2. 2.

Breed $x_{n+1},\ldots,x_{n+m}$ by sequential exchangeable sampling; that is,
let $x_{n+i}\sim P_{\xi}(\cdot\mid x_{1},\ldots,x_{n+i-1})$ for $1\leq i\leq
m$

  3. 3.

Assign ‘survival tickets’ to each of $x_{1},\ldots,x_{n}$; the newly bred
genomes $x_{n+1},\ldots,x_{n+m}$ have as yet no survival tickets.

  4. 4.

Repeat $t$ times:

    1. (a)

Uniformly sample from the population a genome $x_{i}$ which currently has a
survival ticket, and a genome $x_{j}$ which currently does not have a survival
ticket.

    2. (b)

Hold a tournament between $x_{i}$ and $x_{j}$; $x_{i}$ wins with probability
$\frac{w(x_{i})}{w(x_{i})+w(x_{j})}$.

    3. (c)

The winner of the tournament gets the survival ticket; after the tournament,
the winner has the ticket and the loser does not.

  5. 5.

After $t$ tournaments have taken place, the $n$ genomes currently holding
survival tickets are selected to be the new population
$x^{\prime}_{1},\ldots,x^{\prime}_{n}$; the genomes that are not holding
tickets are discarded.

The arguments of section 3.1 can readily be extended to show that this
algorithm has the same stationary distribution (1.2).

###  3.4 Invariance to fitness noise

In real life, survival depends not just on fitness but also on luck. Suppose
that when each new ‘genome’ is ‘born’, it combines its own intrinsic fitness
with an independent random amount of luck. It keeps this amount of luck,
unchanged, throughout its life. Does individual luck at birth alter the
stationary distribution?

To formalise the question, let $\psi_{1},\psi_{2},\ldots,\psi_{n}$ be discrete
random variables that represent ‘luck’. We consider the $\psi_{i}$ as discrete
random variables because the formal derivations become far simpler. The random
variable $\psi_{i}$ can be interpreted as the individual luck that multiplies
the fitness of $i$-th individual. The log fitness function of $i$-th
individual is $\phi(x)+\psi_{i}$.

$P_{n}(\mathbf{x}|\psi_{1},\ldots,\psi_{n})=\frac{1}{Z^{\prime}_{n}}P_{\xi}(\mathbf{x})\exp\big{[}-\sum_{i=1}^{n}\big{(}\phi(x_{i})+\psi_{i}\big{)}\Big{]}$

The joint probability

$\displaystyle P_{n}(\mathbf{x};\psi_{1},\ldots,\psi_{n})$
$\displaystyle=\frac{1}{Z^{\prime}_{n}}P_{\xi}(\mathbf{x})\exp\big{[}-\sum_{i=1}^{n}\big{(}\phi(x_{i})+\psi_{i}\big{)}\big{]}P(\psi_{1},\ldots,\psi_{n})$
$\displaystyle=\frac{1}{Z^{\prime}_{n}}P_{\xi}(\mathbf{x})\exp\big{[}-\sum_{i=1}^{n}\phi(x_{i})\big{]}\exp\big{[}-\sum_{i}\psi_{i}\big{]}P(\psi_{1},\ldots,\psi_{n}).$

We see that the sum factorizes, and so $Z^{\prime}_{n}=Z_{n}\cdot
Z_{n}^{\psi}$, where

$Z_{n}=\sum_{\mathbf{x}}P_{\xi}(\mathbf{x})\exp\big{[}-\sum_{i=1}^{n}\phi(x_{i})\big{]},\quad
Z^{\psi}_{n}=\sum_{\psi_{1},\ldots,\psi_{n}}\exp\big{[}-\sum_{i}\psi_{i}\big{]}\prod_{i}P(\psi_{1},\ldots,\psi_{n})$

Thus the joint measure factorizes

$P_{n}(\mathbf{x};\psi_{1},\ldots,\psi_{n})=P_{n}(\mathbf{x})P_{n}^{\psi}(\psi_{1},\ldots,\psi_{n}),$

where

$P_{n}^{\psi}(\psi_{1},\ldots,\psi_{n})={1\over
Z_{n}^{\psi}}P(\psi_{1},\ldots,\psi_{n})\exp[-\sum_{i}\psi_{i}]$

is a probability measure. Thus, after summing $\psi_{1},\ldots,\psi_{n}$ out,
we end with $P_{n}(\mathbf{x})$:

$\sum_{\psi_{1},\ldots,\psi_{n}}P_{n}(\mathbf{x};\psi_{1},\ldots,\psi_{n})=P_{n}(\mathbf{x}).$

It follows that the stationary distribution is unchanged by multiplicative
noise that is independent of the genomes, for all population sizes. Also,
considering $\psi_{1},\ldots,\psi_{n}$ as the prior, we see that posterior
measure is $P_{n}(\psi_{1},\ldots,\psi_{n})$ that is independent of
$\mathbf{x}$.

###  3.5 Fitness as likelihood: a connection with non-parametric Bayesian MCMC

In Dirichlet Process mixture models, as described for example in [14, 16],
items of data $d_{1},\ldots,d_{n}$ are given, together with a likelihood
function $l(x,d)=P(d\mid x)$, where $x\in\mathcal{X}$ is a discrete latent
variable. The prior distribution over latent variables is given by the
exchangeable process $\xi$, each item of data is associated with its
corresponding latent variable, and the aim of MCMC fitting is to sample latent
variables $x_{1},\ldots,x_{n}$ from the distribution

$P_{n}(x_{1},\ldots,x_{n})=\frac{1}{Z_{n}}P_{\xi}(x_{1},\ldots,x_{n})P(d_{1}\mid
x_{1})\cdots P(d_{n}\mid x_{n})$ (3.4)

where $Z_{n}$ is an appropriate normalizing constant. If we write

${w_{j}(x_{i}):=P(d_{j}\mid x_{i})}$

then this distribution (3.4) is produced by the MCMC algorithm of section 3.1
with the slight modification that the tournament between the new ‘genome’
$x_{n+1}$ and the randomly selected $x_{j}$ is a victory for $x_{n+1}$ with
probability $\frac{w_{j}(x_{n+1})}{w_{j}(x_{j})+w_{j}(x_{n+1})}$.

One might construct an evolutionary “Just-So” story as follows. On a rock in
the ocean, there are $n$ niches, in each of which one member of the species
$\xi$ can live; the fitness of $x$ living in niche $j$ is $w_{j}(x)$.
Evolution occurs when a new individual, bred by ‘n-way recombination’ of all
$n$ parents, then challenges the occupant of a randomly chosen niche by
fighting a tournament, after which the victor survives and takes over the
niche. This gives a (fanciful) evolutionary interpretation to Bayesian latent
variable models with exchangeable priors.

###  3.6 Differences from previous models

The well known Moran Process was introduced by [12], but since then the term
has broadened to include many other related genetic models. Our model has
significant differences from the Moran process and its subsequent variants,
which we now describe.

First, our model implicitly includes mutation; the Moran model did not and the
mutation had to be incorporated independently. To show the difference more
clearly, consider the Polya urn scheme with an $\alpha_{i}>0$ “prior” of
initial balls of color $i$. Unlike in the Moran model, in our model the urn is
not the same as the population. A ball is randomly chosen and returned to the
urn with another ball of the same color. This procedure is repeated $n$ times,
and after that there are $n+|\alpha|$ balls in the urn, but $n$ elements in
the population. Now, according to our breeding rule, a random ball is chosen
again and returned with another ball of the same color. Observe that the prior
balls are involved in the breeding process as much as the rest of the balls.
Whenever a prior ball is chosen during the breeding process, we call it a
mutation. For example, if all $n$ population balls are currently white, and a
black prior ball is chosen, then a black ball enters the population. Or it
might be that a color (type) that has never observed before (inside the $n$
balls taken out) suddenly appears. Or a mutation may produce a colour that is
already common in the population. The main difference from the Moran model is
that the prior balls do not take part in the selection step. Thus, before
tournament(s) start(s), $|\alpha|$ balls ( $\alpha_{i}$ balls from $i$-th
color) are removed from the urn, and so they will never be discarded. There
are now $n+1$ balls in the urn and in the population, one of them will be
ejected according to our selection rule. Now the breeding starts, but before
that the previously removed $|\alpha|$ prior balls are put back to the urn. In
this way, the Markov chain does not have absorbing states, and no fixation
occurs.

Second, in the Moran model and in variants we know of, differences in fitness
are differences in the probability of being selected to breed, whereas in our
model fitness is related to the probability of being selected for discard.
This decision avoids complicating the model of breeding by including arbitrary
fitness, which simplifies the analysis.

Third, our model applies to any exchangeable distributions over finite sets of
possible genomes. The most important instantiation we give here is the product
of Dirichlet processes with finite support, but other complex finitely
supported exchangeable distributions are possible. This means that our model
can be applied to non-trivial genetic algorithms. The Moran model, in
contrast, concerns the fixation probabilities of individual alleles.

Finally, the Moran model required one individual to be born and one to die in
each generation, whereas our model also to applies to exchangeable breeding of
arbitrary numbers of offspring, as described in section 3.3 above.

Our model is a population generalisation of Metropolis Hastings, and is more
closely related to nonparametric Bayes inference with distributions based on
Dirichlet processes, than it is to the Moran model.

Forms of reversible genetic algorithm, or ‘evolutionary MCMC’ that satisfy
detailed balance were previously suggested by [15], [17], and others, but in
these models the genetic operators are not biologically relevant and although
they are described as ‘evolutionary’, they have no relevance to modelling
biological evolution; rather, they are proposal heuristics for Metropolis
Hastings. These methods are normally applied to continuous domains, but if
they were applied to a discrete vector space, the stationary distribution in
our notation would be:

$P_{n}(x_{1},\ldots,x_{n})=\frac{1}{Z_{n}}w(x_{1})\cdots w(x_{n})$

which differs from our stationary distribution in equation 1.2 in that the
breeding term is absent or equivalently that $P_{\xi}$ is the uniform
distribution over $\mathcal{X}$. In our examples, and in genetic models,
$P_{\xi}$ is very far from the uniform distribution.

###  3.7 Designing a reversible evolutionary model

Our larger aim is to develop a model of evolution that is sufficiently
realistic to capture some of the computational power of natural evolution, but
which is also simple and tractable for analysis. To ensure that our model
satisfies detailed balance and has a stationary distribution that factorises
into a breeding and a selection term, we have made the following simplifying
assumptions in addition to the usual simplifications of population genetics or
genetic algorithms:

Overlapping populations

: We believe that overlapping populations are necessary for reversibility. If
full replacement of the population is enforced at each generation, there can
be do guarantee that the population at time $t$ could be easily bred from the
population at time $t+1$. In MCMC, state changes typically occur through
proposing changes that may or may not be accepted; in an overlapping
generations model, if a proposed change is not accepted, we continue with the
same population as before.

$n$-way recombination

: in the product of Dirichlet processes breeding system, each new genome is
bred from $n$ parents rather than from two parents selected from the
population. We conjecture that this is necessary for exact reversibility
because, with long genomes and many mutations, if a child is bred from two
parents, then the child will be more similar to each of its parents than to
other individuals in the population, so that ‘triples’ of two parents and one
child will be identifiable even in the stationary distribution. This breaks
reversibility since the direction of time can be determined observing
evolution in the stationary distribution.

Mutation as sampling

: We consider mutation as sampling from a base distribution of possible
alleles. This model of mutation is not as general as those found in biology,
where mutation probabilities are not symmetric or reversible.

Fitness as lifetime

: All members of the population ‘breed’ at the same rate, and differences in
fitness affect only the expected lifetime of an individual. This clearly
differs from many types of natural selection, but it is also well known that
many organisms continue to produce offspring throughout their lives, so that
their total reproductive success depends on their lifetime, as well as on
other factors.

It is beyond the scope of this article to argue further whether our model
successfully abstracts some essential computational aspects of evolution with
sexual reproduction: we present it merely as an abstraction of sexual
evolution which is significantly more tractable to analyse than other
apparently simple models.

##  4 The measure $P_{n}$

The measure $P_{n}$ is our main object of interest. In this section we show
that the marginal distributions $P_{n}$ over the set of genotypes converge as
the population size $n\rightarrow\infty$; in the next section we characterize
these limits.

Since $\xi$ is exchangeable, by de Finetti’s theorem there exists a prior
measure $\pi$ on the set of all probability measures on ${\cal X}$ (simplex)
${\cal P}$ such that for every $\mathbf{x}\in{\cal X}^{n}$

$P_{\xi}(\mathbf{x})=\int_{{\cal P}}\prod_{i=1}^{n}q(x_{i})\pi(dq)$ (4.1)

so that we can write $P_{n}$ as follows

$P_{n}(\mathbf{x})={1\over Z_{n}}\prod_{i=1}^{n}w(x_{i})\int_{{\cal
P}}\prod_{i=1}^{n}q(x_{i})\pi(dq)={1\over Z_{n}}\int_{{\cal
P}}\prod_{i=1}^{n}\big{(}q(x_{i})w(x_{i})\big{)}\pi(dq),$ (4.2)

where $Z_{n}$ is the normalizing constant. In order to analyze the measure, it
is convenient to rewrite it as follows. First, let us introduce some notation

$\langle q,w\rangle:=\sum_{k=1}^{K}q(k)w(k),\quad
r_{q}(k):={w(k)q(k)\over\langle q,w\rangle},\quad k=1,\ldots,K.$

Note that since $w(k)>0$ for all $k$, $r_{q}$ is correctly defined for every
$q\in\mathcal{P}$. Thus $\langle q,w\rangle$ is the expected weight (under
$q$-measure) and $r_{q}$ is a probability measure on ${\cal P}$. Now

$P_{n}(\mathbf{x})={1\over Z_{n}}\int_{{\cal
P}}\prod_{i=1}^{n}\big{(}q(x_{i})w(x_{i})\big{)}\pi(dq)={1\over
Z_{n}}\int_{{\cal P}}\langle
q,w\rangle^{n}\prod_{i=1}^{n}r_{q}(x_{i})\pi(dq).$ (4.3)

Since $r_{q}$ is a probability measure, it holds that
$\sum_{\mathbf{x}}\prod_{i=1}^{n}r_{q}(x_{i})=1$ and so the normalization
constant for $P_{n}$ is

$Z_{n}=\sum_{\mathbf{x}\in{\cal X}^{n}}\int_{{\cal P}}\langle
q,w\rangle^{n}\prod_{i=1}^{n}r_{q}(x_{i})\pi(dq)=\int_{{\cal P}}\langle
q,w\rangle^{n}\pi(dq).$

Finally note that (4.3) can be rewritten more neatly by defining the measure

$d\bar{\pi}_{n}:={\langle q,w\rangle^{n}\over Z_{n}}d\pi.$ (4.4)

With $\bar{\pi}_{n}$, we have

$P_{n}(\mathbf{x})=\int_{{\cal
P}}\prod_{i=1}^{n}r_{q}(x_{i})\bar{\pi}_{n}(dq)\quad\forall\quad\mathbf{x}\in{\cal
X}^{n}.$ (4.5)

From(4.5), it is easy to find all marginal distributions, namely for any
$m=1,\ldots,n$

$P_{n}(x_{1},\ldots,x_{m})=\int_{{\cal
P}}\prod_{i=1}^{m}r_{q}(x_{i})\bar{\pi}_{n}(dq)={1\over Z_{n}}\int_{{\cal
P}}\langle q,w\rangle^{n-m}\prod_{i=1}^{m}q(x_{i})w(x_{i})\pi(dq).$ (4.6)

In particular, when $(X_{1},\ldots,X_{n})\sim P_{n}$, then

$\displaystyle P(X_{i}=k)$ $\displaystyle=\int_{{\cal
P}}r_{q}(k)\bar{\pi}_{n}(dq)={1\over Z_{n}}\int\langle
q,w\rangle^{n-1}w(k)q(k)\pi(dq)$ $\displaystyle P(X_{i}=k,X_{j}=l)$
$\displaystyle=\int_{{\cal P}}r_{q}(k)r_{q}(l)\bar{\pi}_{n}(dq)={1\over
Z_{n}}\int\langle q,w\rangle^{n-2}w(k)q(k)w(l)q(l)\pi(dq)$

and so on. It is important to observe that $P_{n}(x_{1},\ldots,x_{m})$ depends
on $n$.

##### The limit process.

We have defined for every $n$ the measure (4.5) that describes the genotype
distribution of a $n$-element population. Now the natural question is: do
these measures converge (in some sense) if the population size $n$ grows?
First we have to define the sense of convergence. Since every measure $P_{n}$
is defined on different domain (${\cal X}^{n}$), we cannot speak about
standard (weak) convergence of measures. Instead, we ask about the existence
of a limiting stochastic process. To explain the sense of convergence,
consider that we have defined a triangular array of random variables:

$\displaystyle X_{1,1}\sim P_{1}$ $\displaystyle(X_{2,1},X_{2,2})\sim P_{2}$
$\displaystyle(X_{3,1},X_{3,2},X_{3,3})\sim P_{3}$ $\displaystyle\cdots$
$\displaystyle(X_{n,1},X_{n,2},\ldots,X_{n,n})\sim P_{n}$
$\displaystyle\cdots$

We also know that the joint distribution of the first $m$ variables in every
row depends on $n$. Therefore we ask: is there a stochastic process
$X_{1},X_{2},\ldots$ so that for every $m$ the following convergence holds

$(X_{1,n},\ldots X_{m,n})\Rightarrow(X_{1},\ldots X_{m})?$ (4.7)

According to Kolmogorov’s existence theorem, the existence of a stochastic
process is equivalent to the existence of (finite dimensional) measures
$P^{*}_{m}$ on set ${\cal X}^{m}$, $m=1,2,\ldots$ that satisfy the following
consistency conditions: for every $m$ and for every
$(x_{1},\ldots,x_{m})\in{\cal X}^{m}$, it holds that

$\sum_{x_{m+1}}P_{m+1}^{*}(x_{1},\ldots,x_{m},x_{m+1})=P_{m}^{*}(x_{1},\ldots,x_{m}).$

If we also want (4.7) to be true, then for every $m$ and for every
$(x_{1},\ldots,x_{m})\in{\cal X}^{m}$ the following convergences must hold:

$P_{n}(x_{1},\ldots,x_{m})\to P^{*}(x_{1},\ldots,x_{m}),\quad\forall
m,\quad\forall(x_{1},\ldots,x_{m})\in{\cal X}^{m}.$ (4.8)

We now present a general lemma that guarantees the convergence (4.8). To
achieve the full generality, we let $w$ also depend on $n$. Thus, we have
weights $w_{n}$, and we define the measures $r_{q,n}$ as follows

$r_{q,n}(k):={w_{n}(k)q(k)\over\langle q,w_{n}\rangle}\quad\forall k\in{\cal
X},.$

We start with the following observation, proven in appendix.

###### Claim 4.1

If $w_{n}(i)\to w(i)$ $\forall i\in{\cal X}$, and $r_{q,n}$ and $r_{q}$ are
defined with respect of $w_{n}$ and $w$, respectively, then the following
uniform convergence holds.

$\sup_{q\in{\cal P}}|r_{q,n}(k)-r_{q}(k)|\to 0.$ (4.9)

In the following lemma, $\bar{\pi}_{n}$ is an arbitrary probability measures
on ${\cal P}$, not necessarily as in (4.4). The measure $\bar{\pi}_{n}$ define
$P_{n}$ as in (4.5).

###### Lemma 4.1

Let $w_{n}(k)\to w(k)$ for every $k\in{\cal X}$. If there exists an
probability measure $\bar{\pi}$ such that $\bar{\pi}_{n}\Rightarrow\bar{\pi}$,
then for every $m$ there exists a probability measure $P_{m}^{*}$ on ${\cal
X}^{m}$ so that (4.8) holds. Moreover, for every $(x_{1},\ldots,x_{m})\in{\cal
X}^{m}$

$P^{*}_{m}(x_{1},\ldots,x_{m})=\int_{{\cal
P}}\prod_{i=1}^{m}r_{q}(x_{i})\bar{\pi}(dq),\quad{\rm where}\quad
r_{q}(k)={w(k)q(k)\over\langle w,q\rangle},\quad k=1,\ldots,K$

and the measures $P^{*}_{m}$, $m=1,2,\ldots$ satisfy consistency conditions.

Proof. For every $x_{1},\ldots,x_{m}$ from (4.9), it follows that

$\sup_{q\in{\cal
P}}|\prod_{i=1}^{m}r_{q,n}(x_{i})-\prod_{i=1}^{m}r_{q}(x_{i})|\to 0.$ (4.10)

Since the functions

$q\mapsto\prod_{i=1}^{m}r_{q,n}(x_{i}),\quad
q\mapsto\prod_{i=1}^{m}r_{q}(x_{i})$

are bounded (by 1) continuous functions, from uniform convergence, it holds

$P_{n}(x_{1},\ldots,x_{m})=\int_{{\cal
P}}\prod_{i=1}^{m}r_{q,n}(x_{i})\bar{\pi}_{n}(dq)\to\int_{{\cal
P}}\prod_{i=1}^{m}r_{q}(x_{i})\bar{\pi}(dq)=P^{*}_{m}(x_{1},\ldots,x_{m}).$

Clearly $P_{m}^{*}$ are probability measures. The consistency condition
trivially holds, because

$\sum_{m+1}P^{*}_{m+1}(x_{1},\ldots,x_{m+1})=\int_{{\cal
P}}\sum_{x_{m+1}}\prod_{i=1}^{m+1}r_{q}(x_{i})\bar{\pi}(dq)=\int_{{\cal
P}}\prod_{i=1}^{m}r_{q}(x_{i})\bar{\pi}(dq)=P^{*}_{m}(x_{1},\ldots,x_{m}).$

###  4.1 Frequencies: the measure $Q_{n}$

We now consider how to express the limit measure $P^{*}$ in terms of a
limiting measure on the simplex $\mathcal{P}$. Recall $n_{k}(\mathbf{x})$
defined in (2.1) and let

$\mathbf{n}(\mathbf{x}):=(n_{1},\ldots,n_{K}),\quad\text{where
$n_{i}=n_{i}(\mathbf{x})$, for $i=1,\ldots,K$}$

Since $\xi$ is exchangeable, the probability $P_{\xi}(\mathbf{x})$ depends on
the counts $\mathbf{n}(\mathbf{x})$ only, so

$P_{\xi}(\mathbf{x})=:g(n_{1},\cdots,n_{K}).$

We may now write:

$g(n_{1},\dots,n_{K})=P_{\xi}(\mathbf{x})=\int_{{\cal
P}}\prod_{i=1}^{n}q(x_{i})\pi(dq)=\int_{{\cal
P}}\prod_{k=1}^{K}q(k)^{n_{k}}\pi(dq).$ (4.11)

In what follows, let us denote

$\mathbb{N}_{n}:=\\{(n_{1},\ldots,n_{K}):\sum_{i}n_{i}=n\\}.$

Observe that

$\sum_{(n_{1},\ldots,n_{k})\in\mathbb{N}_{n}}{n!\over n_{1}!\cdots
n_{K}!}g(n_{1},\cdots,n_{K})=1.$

Therefore, the measure $P_{n}$ can be defined on the set $\mathbb{N}_{n}$ as
follows:

$\displaystyle P_{n}({\bf n})$ $\displaystyle:={1\over
Z_{n}}\frac{n!}{n_{1}!\cdots
n_{K}!}g(n_{1},\ldots,n_{K})\prod_{k=1}^{K}\big{(}w_{n}(k)\big{)}^{n_{k}}={n!\over
n_{1}!\cdots n_{K}!}\int_{{\cal
P}}\prod_{k=1}^{K}\big{(}r_{q,n}(k)\big{)}^{n_{k}}\bar{\pi}_{n}(dq).$ (4.12)

Considering the frequencies instead of counts, we can define the corresponding
measure on the simplex ${\cal P}$. Let us denote that measure as $Q_{n}$, so
that with ${\bf n}/n:=(n_{1}/n,\ldots n_{k}/n)$

$Q_{n}\biggl{(}{{\bf n}\over n}\biggr{)}:=P_{n}({\bf n}),\quad\forall{\bf
n}\in\mathbb{N}_{n}.$ (4.13)

Thus $Q_{n}$ is a discrete measure

$Q_{n}=\sum_{{\bf n}\in\mathbb{N}_{n}}P_{n}({\bf n})\delta_{{\bf n}\over n}.$

The advantage of $Q_{n}$ over the measure $P_{n}$ on $\mathbb{N}_{n}$ is that
for any $n$, $Q_{n}$ is defined on the same domain ${\cal P}$, and so one can
speak about the weak convergence of $Q_{n}$. Essentially, obviously, the
measure $P_{n}$ on ${\cal X}^{n}$, the measure $P_{n}$ on $\mathbb{N}_{n}$ and
$Q_{n}$ on ${\cal P}$ are all the same, just the domains are different.

Since the measures $Q_{n}$ are defined on the same space (simplex), it is now
natural to ask, whether there exists a probability measure $Q^{*}$ so that
$Q_{n}\Rightarrow Q^{*}$? It turns out the if the assumption of Lemma 4.1
holds, i.e. $\pi_{n}\Rightarrow{\bar{\pi}}$ and $w_{n}\to w$ (pointwise), then
the limit measure is actually ${\bar{\pi}}r^{-1}$, where

$r:{\cal P}\mapsto{\cal P},\quad r(q)=r_{q}$

and $r$ is defined with respect to limit weight function $w$. Thus for a
measurable $E\subset{\cal P}$,

${\bar{\pi}}r^{-1}(E)={\bar{\pi}}\big{(}r^{-1}(E)\big{)}.$

For example, if ${\bar{\pi}}=\delta_{q^{*}}$ (the measure is concentrated on
one point), then

${\bar{\pi}}r^{-1}=\delta_{r(q^{*})},$

because

$\delta_{q^{*}}r^{-1}(E)=1\quad\Leftrightarrow\quad q^{*}\in
r^{-1}(E)\quad\Leftrightarrow\quad r(q^{*})\in E.$

The following lemma is the counterpart of Lemma 4.1. Again, $\bar{\pi}_{n}$ is
an arbitrary sequence of probability measures on ${\cal P}$, $P_{n}$ are
defined via $\bar{\pi}_{n}$ by (4.12) and $Q_{n}$ via $P_{n}$ as in (4.13).

###### Lemma 4.2

Let $w_{n}(k)\to w(k)$ for every $k\in{\cal X}$. If there exists a probability
measure $\bar{\pi}$ such that $\bar{\pi}_{n}\Rightarrow\bar{\pi}$, then
$Q_{n}\Rightarrow{\bar{\pi}}r^{-1}$.

Proof. Let $f:{\cal P}\to\mathbb{R}$ be a $(K$-variable) bounded continuous
function. By definition of the weak convergence, it suffices to show that

$\int f(q)Q_{n}(dq)\to\int f(q){\bar{\pi}}r^{-1}(dq)=\int
f(r_{q}){\bar{\pi}}(dq),$ (4.14)

where the last equality holds by the change of variable formula. Note that

$\displaystyle\int f(q)Q_{n}(dq)$ $\displaystyle=\sum_{{\bf
n}\in\mathbb{N}_{n}}f\Big{(}{{\bf n}\over n}\Big{)}P_{n}({\bf n})=\sum_{{\bf
n}\in\mathbb{N}_{n}}f\Big{(}{{\bf n}\over n}\Big{)}{n!\over n_{1}!\cdots
n_{K}!}\int\prod_{k=1}^{K}\big{(}r_{q,n}(k)\big{)}^{n_{k}}\bar{\pi}_{n}(dq)$
$\displaystyle=\int\Big{(}\sum_{{\bf n}\in\mathbb{N}_{n}}f\Big{(}{{\bf n}\over
n}\Big{)}{n!\over n_{1}!\cdots
n_{K}!}\prod_{k=1}^{K}\big{(}r_{q,n}(k)\big{)}^{n_{k}}\Big{)}\bar{\pi}_{n}(dq)$
$\displaystyle=\int f_{n}(r_{q,n}(1),\ldots,r_{q,n}(K))\bar{\pi}_{n}(dq)=\int
f_{n}(r_{q,n})\bar{\pi}_{n}(dq),$

where

$f_{n}(r_{q,n}):=:f_{n}(r_{q,n}(1),\ldots,r_{q,n}(K)):=\sum_{(n_{1},\ldots,n_{K})\in\mathbb{N}_{n}:}f\Big{(}{n_{1}\over
n},\ldots,{n_{1}\over n}\Big{)}{n!\over n_{1}!\cdots
n_{K}!}\prod_{k=1}^{K}\big{(}r_{q,n}(k)\big{)}^{n_{k}}$

is the Bernstein polynomial evaluated at
$r_{q,n}=(r_{q,n}(1),\ldots,r_{q,n}(K))$. It is easy to see and well known
that for any vector $r\in{\cal P}$, $f_{n}(r)\to f(r)$, moreover, the
convergence is uniform over ${\cal P}$:

$\sup_{r\in{\cal P}}\Big{|}f_{n}(r)-f(r)\Big{|}\to 0.$ (4.15)

Since for every $q$, $r_{q,n}$ is a probability vector, then

${n!\over n_{1}!\cdots
n_{K}!}\prod_{k=1}^{K}\big{(}r_{q,n}(k)\big{)}^{n_{k}}\leq 1,$

and since $f$ is bounded, we see that for every $n$, $q\mapsto
f_{n}(r_{q,n})=:b_{n}(q)$ is a bounded continuous function. Also the function
$q\mapsto f(r_{q})=:b(q)$ is a bounded continuous function. Then

$\displaystyle\sup_{q}|b_{n}(q)-b(q)|$
$\displaystyle\leq\sup_{q}|f_{n}(r_{q,n})-f(r_{q,n})|+\sup_{q}|f(r_{q,n})-f(r_{q})|.$

By (4.9), $\sup_{q}|r_{q,n}(i)-r_{q}(i)|\to 0$ for every $i$. Then also
$\sup_{q}\|r_{q,n}-r_{q}\|\to 0$. A continuous function on compact space is
uniformly continuous, so

$\sup_{q}|f(r_{q,n})-f(r_{q})|\to 0.$

By (4.15),

$\sup_{q}|f_{n}(r_{q,n})-f(r_{q,n})|\leq\sup_{r}|f_{n}(r)-f(r)|\to 0.$

Therefore, we have shown that $\sup_{q}|b_{n}(q)-b(q)|$ implying that

$\int f(q)Q_{n}(dq)=\int b_{n}(q)\bar{\pi}_{n}(dq)\to\int
b(q){\bar{\pi}}(dq)=\int f(r_{q}){\bar{\pi}}(dq).$

##  5 $P^{*}$ and $Q^{*}$ in the large population limit

In what follows, let us rewrite fitnesses in terms of $\phi(k):=-\ln(w(k))$,
so that for any genome $k\in{\cal X}$, $w(k)=\exp(-\phi(k))$. Moreover, in
order to increase the influence of prior, we let weights $w_{n}$ depend on $n$
in the following way:

$w_{n}(k)=\exp[-{\phi(k)\over n^{\lambda}}],\quad k=1,\ldots,K$ (5.1)

where $\lambda\geq 0$ and $0\leq\phi(1)<\phi(2)<\cdots<\phi(K)$. The case
$\lambda=0$ corresponds to fitness that is constant in that it does not vary
with $n$. Clearly, for every $k$, $w_{n}(k)\to w(k)$ and $\lambda$ controls
the speed of that convergence. When $\lambda>0$, then $w(k)=1$ implying that
in this case the mapping $r$ is identity, i.e. for every $q$, $r_{q}=q$.

Let us return to our original $\bar{\pi}_{n}$, defined as in (4.4) with
$w_{n}$:

$\bar{\pi}_{n}(E)=\int_{E}{\langle q,w_{n}\rangle^{n}\over Z_{n}}\pi(dq),\quad
Z_{n}=\int\langle q,w\rangle^{n}\pi(dq).$ (5.2)

In this section we consider the case where the prior measure $\pi$ is
independent of $n$, and the support of $\pi$ is the whole simplex ${\cal P}$.
Since by assumption $w(1)>w(2)$, clearly the function $q\mapsto\langle
q,w\rangle$ has unique maximizer $q^{*}:=(1,0,\ldots,0)$. The following
theorem states that phase transition occurs: when $0\leq\lambda<1$, then
$\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$ and $P_{n}(x_{1},\ldots,x_{m})\to
P^{*}(x_{1},\ldots,x_{m}),$ where

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}r_{q^{*}}(x_{i})=\prod_{i=1}^{m}q^{*}(x_{i})=\left\\{\begin{array}[]{ll}1,&\hbox{if
for every $i$, $x_{i}=1$ ;}\\\ 0,&\hbox{else.}\end{array}\right.,$

because in both cases (i.e. $\lambda=0$ and $\lambda\in(0,1)$), it holds that
$r_{q^{*}}=q^{*}$. Thus the limit process $X_{1},X_{2},\ldots$ has only one
realization: $1,1,\ldots$. In this case also $Q_{n}\Rightarrow\delta_{q^{*}}$.
Thus, when $\lambda\in[0,1)$ then only the fittest genotype survives, no one
else has any change, no matter what the prior says. In other words, the
influence of the prior vanishes.
When $\lambda=1$, then $\bar{\pi}_{n}$ as well $Q_{n}$ converges to a
nondegenerate distribution, specified below, and also the limit measure
$P^{*}$ is non-degenerate.
And finally, when $\lambda>1$, then $\bar{\pi}_{n}\Rightarrow\pi$,
$Q_{n}\Rightarrow\pi$ and the measure $P^{*}$ is the law of birth process
$\xi$. In this case the influence of fitness vanish and only the prior matters
– the limit process equals to the breeding one.

###### Theorem 5.1

Let the fitness function be defined as in (5.1) and assume that the support of
the prior $\pi$ is ${\cal P}$. Then the following convergences hold:

1)

If $\lambda\in[0,1)$, then $\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$,
$Q_{n}\Rightarrow\delta_{q^{*}}$ and (4.8) holds with

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}q^{*}(x_{i}),\text{ where
}q^{*}=(1,0,\ldots,0).$

2)

If $\lambda=1$, then $\bar{\pi}_{n}\Rightarrow{\bar{\pi}}$,
$Q_{n}\Rightarrow{\bar{\pi}}$ and (4.8) holds with

$P^{*}(x_{1},\ldots,x_{m})=\int\prod_{i=1}^{m}q(x_{i})\bar{\pi}(dq),$

where for every $E\subset{\cal P}$,

${\bar{\pi}}(E)={1\over Z}\int_{E}\exp[-\langle\phi,q\rangle]\pi(dq),\quad
Z=\int\exp[-\langle\phi,q\rangle]\pi(dq).$

3)

If $\lambda>1$, then $\bar{\pi}_{n}\Rightarrow{\pi}$, $Q_{n}\Rightarrow{\pi}$
and (4.8) holds with

$P^{*}(x_{1},\ldots,x_{m})=P_{\xi}(x_{1},\ldots,x_{m}).$

###  5.1 Proof of Theorem 5.1

Before proving the theorem, let us state a very useful preliminary result.
Recall that simplex ${\cal P}$ is a compact set. Let $f_{n},f:{\cal
P}\to\mathbb{R}^{+}$ be continuous, hence bounded measurable functions so that
$f_{n}\to f$ uniformly and let $m_{n}\to\infty$ be an increasing sequence. We
are given a measure $\pi$ on ${\cal P}$, and we are interested in the
asymptotic behavior of the measure $\nu_{n}$, where

$\nu_{n}(E):=\int_{E}h_{n}(q)\pi(dq),\quad
h_{n}(q):={f_{n}^{m_{n}}(q)\over\int
f_{n}^{m_{n}}(q)\pi(dq)}=\Big{(}{f_{n}(q)\over\|f_{n}\|_{m_{n}}}\Big{)}^{m_{n}}.$

Here we assume that $\int f_{n}^{m_{n}}d\pi<\infty$ for every $n$. If $\pi$ is
a finite measure, then the conditions automatically holds due to the
boundedness of $f_{n}$. In what follows, let

${\cal S}^{*}:=\\{q\in{\cal S}:f(q)=\|f\|_{\infty}\\},\quad{\cal
S}_{\delta}^{*}:=\\{q\in{\cal S}:f(q)>\|f\|_{\infty}-\delta\\}.$ (5.3)

Here $\|f\|_{\infty}$ is the essential supremum of $f$ with respect to the
$\pi$-measure. If $f$ is continuous and the support of $\pi$ is ${\cal P}$,
then $\|f\|_{\infty}=\sup_{q}f(q)$. The proof of the following Proposition 5.1
is given in appendix.

###### Proposition 5.1

Let $f_{n}\to f$ uniformly and let $\pi$ be a finite measure on ${\cal P}$.
Then for every $\delta>0$, $\nu_{n}\big{(}{\cal S}_{\delta}^{*}\big{)}\to 1$.
If ${\cal S}^{*}=\\{q^{*}\\}$, then $\nu_{n}\Rightarrow\delta_{q^{*}}.$

Besides Proposition 5.1, the proof of Theorem 5.1 is based on the following
well-known observation: when $m\to\infty$, then

$\sup_{q\in{\cal P}}\Big{|}\big{\langle}\exp[-{\phi\over
m}],q\big{\rangle}^{m}-\exp[-\langle\phi,q\rangle]\Big{|}\to 0.$ (5.4)

Proof. (Theorem 5.1)

1)

For $\lambda=0$, take $f_{n}(q)=f(q)=\langle w,q\rangle$. From Proposition
5.1, it follows that $\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$. Since for any
weight $w$, $r_{q^{*}}(k)=q^{*}(k)$, from Lemma 4.1, it follows that

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}r_{q^{*}}(x_{i})=\prod_{i=1}^{m}q^{*}(x_{i}).$

Since $\bar{\pi}_{n}r^{-1}(q^{*})=q^{*}$, from Lemma 4.2, it follows that
$Q_{n}\Rightarrow\delta_{q^{*}}$. So, for $\lambda=0$, the statement is proven
and we now consider the case $\lambda\in(0,1)$. Let

$f_{n}(q):=\langle\exp[-{\phi\over n^{\lambda}}],q\rangle^{n^{\lambda}},\quad
f(q):=\exp[-\langle\phi,q\rangle].$

By (5.4), $\|f_{n}-f\|_{\infty}\to 0$. Since $\lambda\in(0,1)$, take
$m_{n}=n^{1-\lambda}$. Then

$h_{n}(q):={\langle\exp[-{\phi\over n^{\lambda}}],q\rangle^{n}\over Z_{n}}$

is the density of $\bar{\pi}_{n}$ with respect to $\pi$. Since $f$ is
continuous, the set ${\cal S}^{*}$ in (5.3) is

${\cal S}^{*}=\arg\max_{q\in{\cal S}}f(q)=\arg\min_{q\in{\cal
S}}\langle\phi,q\rangle=\\{(1,0,\ldots,0)\\}=\\{q^{*}\\}.$

By Proposition 5.1, $\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$. As in the case
of $\lambda=0$, it follows that $Q_{n}\Rightarrow\delta_{q^{*}}$ and
$P^{*}(x_{1},\ldots,x_{m})=1$ if and only if $x_{1}=\cdots=x_{m}=1$.

2)

Since for any $q$ and any $n$, it holds

$\langle\exp[-{\phi\over n}],q\rangle^{n}\leq e^{-\phi(1)}\leq 1,$

we obtain from (5.4) and bounded convergence that for any measurable $E$

$\int_{E}\langle\exp[-{\phi\over
n}],q\rangle^{n}\pi(dq)\to\int_{E}\exp[-\langle\phi,q\rangle]\pi(dq).$ (5.5)

Recall

$\bar{\pi}(E)={1\over Z}{\int_{E}}\exp[-\langle\phi,q\rangle]\pi(dq),\quad{\rm
where}\quad Z=\int\exp[-\langle\phi,q\rangle]\pi(dq).$

Therefore, from (5.5), it follows that when $\lambda=1$, we have

$\bar{\pi}_{n}(E)\to{\bar{\pi}}(E),$

meaning that $\bar{\pi}_{n}\Rightarrow\bar{\pi}$ (even in a stronger sense).
Since $r_{q}=q$, from Lemma 4.1, it follows that the limits of
$P_{n}(x_{1},\ldots,x_{m})$ are

$P^{*}(x_{1},\ldots,x_{m})=\int\prod_{i=1}^{m}q(x_{i})\bar{\pi}(dq)={1\over
Z}\int\prod_{i=1}^{m}q(x_{i})\exp[-\langle\phi,q\rangle]\pi(dq).$

Since $r$ is identity function, by Lemma 4.2 the limit measure of frequencies
is $\bar{\pi}$, i.e. $Q_{n}\Rightarrow\bar{\pi}$.

3)

Since for any $q$,

$\langle\exp[-{\phi\over n^{\lambda}}],q\rangle^{n}\to 1$

by dominated convergence, again, for any measurable $E$

$\int_{E}\langle\exp[-{\phi\over n}],q\rangle^{n}\pi(dq)\to\pi(E).$

Therefore $\bar{\pi}_{n}\Rightarrow\pi$. By Lemma 4.1, the limits of
$P_{n}(x_{1},\ldots,x_{m})$ are

$P^{*}(x_{1},\ldots,x_{m})=\int\prod_{i=1}^{m}q(x_{i})\pi(dq),$

so that the limit process is $\xi$. The convergence $Q_{n}\Rightarrow\pi$
follows from Lemma 4.2.

We have seen that the critical case $\lambda=1$ is the only case where the
prior and fitnesses both determine the limit measure. In this case, the limit
process $X_{1},X_{2},\ldots$, governed by $P^{*}$ has marginals

$\displaystyle P(X_{i}=k)$ $\displaystyle=P(X_{1}=k)=\int
q(k)\bar{\pi}(dq)={1\over Z}\int e^{-\langle\phi,q\rangle}q(k)\pi(dq),$
$\displaystyle P(X_{i}=k,X_{j}=l)$ $\displaystyle=P(X_{1}=k,X_{2}=l)={1\over
Z}\int e^{-\langle\phi,q\rangle}q(k)q(l)\pi(dq).$

It is also interesting to point out that in the critical case $\lambda=1$, the
measure $\bar{\pi}$ satisfies

$\bar{\pi}=\arg\min_{\pi^{\prime}\in E}D(\pi\|\pi^{\prime}),$

where $E$ is a set of probability measures on ${\cal P}$, namely
$E:=\\{\pi^{\prime}:\int\langle\phi,q\rangle\pi^{\prime}(dq)\geq c\\}$ and
$c>0$ is a constant.

##  6 Dirichlet prior

Also in the current section we consider the weights $w_{n}(k)$ as in (5.1),
where $\lambda\in[0,1]$. We already know that in the case of constant priors,
the case $\lambda<1$ means that the fitnesses will prevail over the prior, and
the limit measure is degenerate one. Therefore, it is meaningful to consider
the non-constant priors so that the influence of prior increases with suitable
rate. Therefore, in the present section, we consider Dirichlet’ priors

$\pi_{n}={\rm Dir}(n^{1-\lambda}\alpha_{1},\ldots n^{1-\lambda}\alpha_{K}),$
(6.1)

where $\alpha:=(\alpha_{1},\ldots,\alpha_{K})$, $\alpha_{k}>0$ and
$|\alpha|:=\sum_{k}\alpha_{k}$. The constant
$\sum_{i}n^{1-\lambda}\alpha_{k}=|\alpha|n^{1-\lambda}$ is the so called
concentration or precision parameter, the bigger that parameter, the more
prior is concentrated over it expectation
$(\alpha_{1}/|\alpha|,\ldots,\alpha_{K}/|\alpha|)$. Increasing the
concentration parameter increases the influence of prior, and now it is clear
the the smaller is $\lambda$, the bigger must be the prior influence. This
justifies the choice of $n^{1-\lambda}$. The case $\lambda=1$ corresponds to
already studied case of constant priors, therefore we now consider the case
$\lambda\in[0,1)$ The following theorem shows that the phase transition occurs
again.

###### Theorem 6.1

Let the fitness function be defined be defined as in (5.1) and the prior
$\pi_{n}$ as in (6.1). Let $\bar{\pi}_{n}$ be defined as in (5.2) with
$\pi_{n}$ instead of $\pi$. Then the following convergences hold:

1)

If $\lambda=0$, then $\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$, where $q^{*}$
is the unique maximizer of the following function

$\ln\langle e^{-\phi},q\rangle+\sum_{k=1}^{K}\alpha_{k}\ln q(k).$ (6.2)

Then $Q_{n}\Rightarrow\delta_{r^{*}}$, where $r^{*}=r_{q^{*}}$, so that
$r^{*}(k)\propto q^{*}(k)w(k)$, and (4.8) holds with

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}r^{*}(x_{i}).$

2)

If $\lambda\in(0,1)$, then $\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$, where
$q^{*}$ is the unique maximizer of the following function:

$-\langle\phi,q\rangle+\sum_{k=1}^{K}\alpha_{k}\ln q(k)$ (6.3)

Then $Q_{n}\Rightarrow\delta_{q^{*}}$ and (4.8) holds with

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}q^{*}(x_{i}).$

Let us start with proving the uniqueness of the solutions of (6.2) and (6.3).
The proof of the following lemma is in appendix.

###### Lemma 6.1

1)

The function (6.2) has an unique maximizer $q^{*}$, where

$q^{*}(k)={\alpha_{k}\over(1+|\alpha|)-{w(k)\over\theta}},\quad k=1,\ldots,K$
(6.4)

where $\theta>0$ is a parameter satisfying $\theta=\langle w,q^{*}\rangle$.

2)

The function (6.3) has an unique maximizer $q^{*}$, where

$q^{*}(k)={\alpha_{k}\over\phi(k)+|\alpha|-\theta},\quad k=1,\ldots,K,$ (6.5)

where $\theta>0$ is the parameter satisfying
$\theta=\langle\phi,q^{*}\rangle.$

##### Proof of theorem 6.1.

1)

In the case $\lambda=0$, the measure $\bar{\pi}_{n}$ has the following density
with respect to the Lebesgue measure:

$\bar{\pi}_{n}(q)={1\over Z_{n}}\langle e^{-\phi},q\rangle^{n}\cdot{1\over
B(n\alpha)}\prod_{k=1}^{K}(q(k))^{n(\alpha_{k}-{1\over n})}={f_{n}(q)^{n}\over
Z^{\prime}_{n}},$

where

$f_{n}(q):=\langle
e^{-\phi},q\rangle\prod_{k=1}^{K}(q(k))^{(\alpha_{k}-{1\over n})},\quad
Z^{\prime}_{n}:=\int f^{n}_{n}(q)dq.$

Clearly for every $q$, $f_{n}(q)\to f(q)$, where

$f(q)=\langle e^{-\phi},q\rangle\prod_{k=1}^{K}(q(k))^{\alpha_{k}}.$

It is not hard to see that the convergence is uniform, i.e.
$\sup_{q}|f_{n}(q)-f(q)|\to 0$. By 1) of Lemma 6.1, the function $f$ has
unique maximizer $q^{*}$ (6.4), i.e. ${\cal S}^{*}=\\{q^{*}\\}$. Now apply
Proposition 5.1 with $\pi$ being the Lebesgue measure on ${\cal P}$ (hence
$\pi$ is finite) and $m_{n}=n$ so that

$\nu_{n}(E)=\int_{E}h_{n}dq={1\over
Z^{\prime}_{n}}\int_{E}{f_{n}(q)^{n}}dq=\bar{\pi}_{n}(E).$

Since all assumptions are fulfilled, we have
$\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$. Since $w_{n}=w$, by Lemma 4.1 the
limit process has finite-dimensional distributions

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}r^{*}(x_{i}),\quad\text{where}\quad
r^{*}=r_{q^{*}}$

so that the limit process $P^{*}$ corresponds to a i.i.d. sequence
$X_{1},X_{2},\ldots$ with $X_{1}\sim r^{*}$. According to Lemma 4.2, the
frequencies $Q_{n}$ converge weakly to the measure ${r^{*}}$ and this is also
quite obvious by SLLN.

2)

The proof is similar: $\bar{\pi}_{n}$ has density (with respect to Lebesgue
measure)

${f_{n}(q)^{m_{n}}\over Z^{\prime}_{n}},\quad\text{where}\quad
f_{n}(q)=\langle e^{-{\phi\over
n^{\lambda}}},q\rangle^{n^{\lambda}}\prod_{k=1}^{K}(q(k))^{(\alpha_{k}-{1\over
n^{1-\lambda}})},\quad m_{n}=n^{(1-\lambda)}.$

Since

$\langle e^{-{\phi\over n^{\lambda}}},q\rangle^{n^{\lambda}}\to
e^{\langle\phi,q\rangle},$

uniformly over $q$, we have that sequence $f_{n}$ converges uniformly to

$f(q)=e^{-\langle\phi,q\rangle}\prod_{k=1}^{K}(q(k))^{\alpha_{k}}.$

By 2) of Lemma 6.1, the function $f$ has unique maximizer $q^{*}$ (6.5). As in
the case 1), it is easy to see that the assumptions of Proposition 5.1 are
fulfilled with $\pi$ being Lebesgue measure on ${\cal P}$,
$m_{n}=n^{1-\lambda}$ and so $\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$. In the
present case, for every $k=1,\ldots,K$, $w_{n}(k)\to 1$ and so by Lemma 4.1,
the limit process $P^{*}$ is i.i.d. process with distribution ${q^{*}}$
(because $r$ is identity function). According to Lemma 4.2, the frequencies
$Q_{n}$ converge weakly to the measure $\delta_{q^{*}}$.

###  6.1 Relation between $\lambda=0$ and $\lambda\in(0,1)$

From (5.4), it follows:

$\sup_{q\in{\cal P}}\Big{|}m\ln\langle\exp[-{\phi\over
m}],q\rangle+\langle\phi,q\rangle]\Big{|}\to 0$ (6.6)

so that with

$f_{m}(q):=\ln\langle\exp[-{\phi\over
m}],q\big{\rangle}+\sum_{k}{\alpha_{k}\over m}\ln q(k),\quad
f(q)=-\langle\phi,q\rangle+\sum_{k}{\alpha_{k}}\ln q(k),$

we have

$\sup_{q\in{\cal P}}|mf_{m}(q)-f(q)|\to 0.$

Since $f(q)$ is as in (6.3), it has the unique maximizer $q^{*}$ given in
(6.4). On the other hand, the maximizer of $mf_{m}(q)$ is the same as the
maximizer of $f_{m}(q)$, which corresponds to (6.2) where $\phi$ is replaced
by $\phi/m$ and $\alpha$ is replaced by $\alpha/m$. Let this unique maximizer
be $q^{*}_{m}$ Since the functions $mf_{m}(\cdot)$ and $f(\cdot)$ are
continuous, uniformly convergent and having unique maximum, it follows that
$q^{*}_{m}\to q^{*}$ (in usual sense, because ${\cal P}$ is compact). Thus, we
have proven the following proposition.

###### Proposition 6.1

Let

$q^{*}_{m}=\arg\max_{q}\Big{(}\ln\langle\exp[-{\phi\over
m}],q\big{\rangle}+\sum_{k}{\alpha_{k}\over m}\ln q(k)\Big{)}$

and let $q^{*}$ be the maximizer of (6.3). Let $r_{m}^{*}$ be the
corresponding $r$ measure, i.e. $r_{m}^{*}(k)\propto
q^{m}(k)\exp[-{\phi(k)\over m}].$ Then $q^{*}_{m}\to q^{*}$ and $r^{*}_{m}\to
q^{*}$.

###  6.2 Product of Dirichlet priors

Recall the setup in Subsection 2.1. The set of genomes is now ${\cal
X}^{L}=\overbrace{{\cal X}\times\cdots\times{\cal X}}^{L}$ and the breeding
process $\xi=(\xi^{1},\ldots,\xi^{L})$, where $\xi^{l}$ are independent
exchangeable processes. We now assume that the prior of $\xi^{l}$ is
$\pi^{l}={\rm Dir}(\alpha^{l}),$ where
$\alpha^{l}=(a_{1}^{l},\ldots,\alpha_{K}^{l})$, $l=1,\ldots,L$. In this model,
$L$ different Polya urns are run independently. Let ${\cal P}^{L}$ be the set
of $L$-fold product measures:

${\cal P}^{L}:=\\{q^{1}\times\cdots\times q^{L}:q^{j}\in{\cal P}\\},$

where ${\cal P}$, as previously, stands for the $(K-1)$-dimensional simplex.
Observe that ${\cal P}^{L}$ is a compact subset of the set of all possible
probability measures on ${\cal X}^{L}$. Since the components of $\xi$ are
independent, the prior $\pi$ of $\xi$ is the product of Dirichlet measures
$\pi=\pi^{1}\times\cdots\times\pi^{L}$. This means that the support of $\pi$
is ${\cal P}^{L}$ and for every element $q=q^{1}\times\cdots\times
q^{L}\in{\cal P}^{L}$, the density is (with slight abuse of notation, $\pi$
stands for the measure as well as for its density)

$\pi(q)=\prod_{l=1}^{L}\pi^{l}(q^{l})={1\over
B}\prod_{l=1}^{L}\prod_{k=1}^{K}\big{(}q(k)^{l}\big{)}^{\alpha^{l}_{k}-1},\quad
B:=\prod_{l=1}^{L}B(\alpha^{l}).$

The function $\phi$ is now defined on the set ${\cal X}^{L}$, and so for any
$q\in{\cal P}^{L}$,

$\langle\phi,q\rangle=\sum_{(k_{1},\ldots,k_{L})\in{\cal
X}^{L}}\phi(k_{1},\ldots,k_{L})q^{1}(k_{1})\cdots q^{L}(k_{L})$

and $\langle e^{-\phi},q\rangle$ is defined similarly. When $\lambda=0$, the
measure $\bar{\pi}_{n}$ has density $f_{n}(q)^{n}/Z^{\prime}_{n}$, where

$f_{n}(q)=\langle
e^{-\phi},q\rangle\prod_{l=1}^{L}\prod_{k=1}^{K}(q^{l}(k))^{(\alpha^{l}_{k}-{1\over
n})},\quad Z^{\prime}_{n}:=\int f^{n}_{n}(q)dq.$

Clearly $f_{n}(q)$ converges uniformly to

$f(q):=\langle
e^{-\phi},q\rangle\prod_{l=1}^{L}\prod_{k=1}^{K}(q^{l}(k))^{(\alpha^{l}_{k})}.$
(6.7)

Similarly, when $\lambda\in(0,1)$ the measure $\bar{\pi}_{n}$ has density
$f_{n}(q)^{m_{n}}/Z^{\prime}_{n}$, where $m_{n}=n^{1-\lambda}$,

$f_{n}(q)=\langle e^{-\phi\over
n^{\lambda}},q\rangle^{n^{\lambda}}\prod_{l=1}^{L}\prod_{k=1}^{K}(q^{l}(k))^{(\alpha^{l}_{k}-{1\over
n^{1-\lambda}})},\quad Z^{\prime}_{n}:=\int f^{m_{n}}_{n}(q)dq.$

Again, $f_{n}(q)$ converges uniformly to

$f(q):=e^{\langle{-\phi},q\rangle}\prod_{l=1}^{L}\prod_{k=1}^{K}(q^{l}(k))^{(\alpha^{l}_{k})}.$
(6.8)

When (6.7) (resp. (6.8)) have unique solution $q^{*}$, then the statements of
Theorem 6.1 hold (the proof is the same):

1)

Suppose $\lambda=0$ and (6.7) has unique maximizer
$q^{*}=q^{1}\times\cdots\times q^{L}$. Then
$\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$, and
$Q_{n}\Rightarrow\delta_{r^{*}}$, where

$r^{*}(k_{1},\ldots,k_{L})\propto w(k_{1},\ldots,k_{L})q^{1}(k_{1})\cdots
q^{L}(k_{L}),$

and $w(k_{1},\ldots,k_{L})=\exp[-\phi(k_{1},\ldots,k_{L})].$ Observe that the
measure $r^{*}$ is not necessarily a product measure. Then also (4.8) holds
with

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}r^{*}(x_{i}),\quad x_{i}\in{\cal
X}^{L}.$

2)

Suppose $\lambda\in(0,1)$ and (6.8) has unique maximizer
$q^{*}=q^{1}\times\cdots\times q^{L}$. Then
$\bar{\pi}_{n}\Rightarrow\delta_{q^{*}}$, $Q_{n}\Rightarrow\delta_{q^{*}}$ and
(4.8) holds with

$P^{*}(x_{1},\ldots,x_{m})=\prod_{i=1}^{m}q^{*}(x_{i}),\quad x_{i}\in{\cal
X}^{L}.$

In the case $L>1$, the maximizer of (6.7) and (6.8) is not always unique.
Whether it is unique or not depends on $\phi$ and vectors $\alpha^{l}$.
Indeed, maximizing (6.7) is equivalent to minimizing

$-\ln\big{(}\langle
e^{-\phi},q\rangle\big{)}-\sum_{l=1}^{L}\sum_{k=1}^{K}\alpha_{k}^{l}\ln(q^{l}(k)),$
(6.9)

and $-\sum_{l=1}^{L}\sum_{k=1}^{K}\alpha_{k}^{l}\ln(q^{l}(k))$ is always a
convex function. So, when the parameters $\alpha^{l}_{k}$ are big enough, then
the whole function (6.9) becomes convex. The same argument holds for (6.8). We
shall present some sufficient conditions for convexity of (6.7) and (6.8) for
the case $K=L=2$ below.

Recall: when a positive continuous function $f(q)$ has one maximizer $q^{*}$,
then for any sequence $m_{n}\to\infty$, the measures $\nu_{n}$ with densities
proportional to $f^{m_{n}}(q)$ converge weakly to $\delta_{q^{*}}$. When the
function has, say, two maximizers, $q^{*}_{1}$ and $q^{*}_{2}$, then by
Proposition 5.1, for all disjoint open balls $B_{1}$ and $B_{2}$ so that
$q^{*}_{i}\in B_{i}$, it still holds that $\nu_{n}(B_{1})+\nu_{n}(B_{2})\to
1$. Thus, when the measures $\nu_{n}$ are weakly convergent, then the limit
measure is concentrated on $\\{q^{*}_{1},q^{*}_{2}\\}$ so that the limit
measure must be $p\delta_{q^{*}_{1}}+(1-p)\delta_{q^{*}_{2}}$, for some
$p\in[0,1]$. In this case $\nu_{n}(B_{1})\to p$. However, the function $f$
might be so that the limits $\nu_{n}(B_{i})$ do not exist. And even if they do
exist (i.e. the measures $\nu_{n}$ are weakly convergent), the limits $p$ and
$1-p$ might be arbitrary real numbers, and hard to determine. Therefore the
following theorem adapted from [9] (Theorem 5.7 and a remark after it) might
be very useful.

###### Theorem 6.2

Suppose $K\subset\mathbb{R}^{k}$ is a compact non-empty subset and let
$g:K\to[0,\infty)$ be a twice continuously differentiable function with
finitely many minimum points $\\{a_{1},\ldots,a_{r}\\}$ all located in the
interior of $K$. Let, for every $i=1,\ldots r$ the Hessian of $g$ at $a_{i}$
be positive definite. Given any increasing sequence $m_{n}$, define the
sequence of measures

$\nu_{n}(E):={1\over Z_{n}}\int_{E}\exp[-m_{n}\cdot g(x)]dx,\quad
Z_{n}:=\int\exp[-m_{n}\cdot g(x)]dx.$

Then $\nu_{n}\Rightarrow\nu$, where

$\nu=\sum_{i=1}^{r}p_{i}\delta_{a_{i}},\quad\text{ with
}p_{i}\propto{1\over\sqrt{\det H(a_{i})}}$

and $\det H(a_{i})$ is a determinant of Hessian evaluated at $a_{i}$.

To apply the theorem in our case, let us first note that any $K-1$ dimensional
simplex can be considered as a $K-1$-dimensional non-empty compact set

${\cal P}_{K}=\\{(q(1),\ldots,q(K-1)):q(k)\geq 0,\sum_{k}^{K-1}q(k)\leq 1\\}.$

Therefore, our search space ${\cal P}^{L}$ can be considered as a subset in
$\mathbb{R}^{L(K-1)}.$ This subset has non-empty interior. Clearly any
solution of (6.7) and (6.8) has all components strictly positive so that all
maximizers of maximizer of (6.7) are interior points and the same holds for
(6.8). We take $g(q)=-\ln f(q)$, where $f$ is as in (6.7) or (6.8). Thus, for
any $m$, $\exp[-mg(y)]=f^{m}(q)$ so that the measure defined in the statement
of theorem is

$\nu_{n}(E)\propto\int_{E}f^{m_{n}}(q)dq.$

However, even when the measures $\nu_{n}$ converge weakly to a limit, it does
not automatically follow that the measures $\bar{\pi}_{n}$ converge to the
same limit even if $f_{n}$ converges to $f$ uniformly. This convergence might
depend on the speed of the uniform convergence, and we leave it for the
further studies and proceed with an example instead.

####  6.2.1 The case $K=L=2$

Let us analyze more closely the case $K=2$ and $L=2$. Denote $q^{1}(1)=:z_{1}$
and $q^{2}(1)=:z_{2}$. Also denote $\alpha_{k}^{1}=\alpha_{k}$ and
$\alpha_{k}^{2}=\beta_{k}$. The function (6.9) is

$g(z_{1},z_{2})=-\ln\big{(}\langle w,z\rangle\big{)}-\alpha_{1}\ln
z_{1}-\alpha_{2}\ln(1-z_{1})-\beta_{1}\ln z_{2}-\beta_{2}\ln(1-z_{2}),$ (6.10)

where

$\langle
w,z\rangle=w(1,1)z_{1}z_{2}+w(1,2)z_{1}(1-z_{2})+w(2,1)(1-z_{1})z_{2}+w(2,2)(1-z_{1})(1-z_{2}).$

Thus with $w^{*}:=w(1,1)-w(1,2)-w(2,1)+w(2,2)$ and

$\displaystyle\theta^{1}_{1}$
$\displaystyle:=w(1,1)z_{2}+w(1,2)(1-z_{2}),\quad\theta^{1}_{2}=w(2,1)z_{2}+w(2,2)(1-z_{2})$
$\displaystyle\theta^{2}_{1}$
$\displaystyle:=w(1,1)z_{1}+w(2,1)(1-z_{1}),\quad\theta^{2}_{2}:=w(1,2)z_{1}+w(2,2)(1-z_{1}).$

we obtain the Hessian

$\displaystyle\left(\begin{array}[]{cc}{\partial^{2}g\over\partial
z_{1}^{2}}&{\partial^{2}g\over\partial z_{1}\partial z_{2}}\\\
{\partial^{2}g\over\partial z_{2}\partial z_{1}}&{\partial^{2}g\over\partial
z_{2}^{2}}\\\ \end{array}\right)=\left(\begin{array}[]{cc}{\alpha_{1}\over
z_{1}^{2}}+{\alpha_{2}\over(1-z_{1})^{2}}+{(\theta_{1}^{1}-\theta_{2}^{1})^{2}\over\langle
w,z\rangle^{2}}&{(\theta_{1}^{1}-\theta_{2}^{1})(\theta_{1}^{2}-\theta_{2}^{2})-w^{*}\langle
w,z\rangle\over\langle w,z\rangle^{2}}\\\
{(\theta_{1}^{1}-\theta_{2}^{1})(\theta_{1}^{2}-\theta_{2}^{2})-w^{*}\langle
w,z\rangle\over\langle w,z\rangle^{2}}&{\beta_{1}\over
z_{2}^{2}}+{\beta_{2}\over(1-z_{2})^{2}}+{(\theta_{1}^{2}-\theta_{2}^{2})^{2}\over\langle
w,z\rangle^{2}}\\\ \end{array}\right).$

The elements in the main diagonal are strictly positive and therefore the
matrix is positive definite if the determinant is positive, i.e.

$\displaystyle\Big{(}{\alpha_{1}\over
z_{1}^{2}}+{\alpha_{2}\over(1-z_{1})^{2}}+{(\theta_{1}^{1}-\theta_{2}^{1})^{2}\over\langle
w,z\rangle^{2}}\Big{)}\Big{(}{\beta_{1}\over
z_{2}^{2}}+{\beta_{2}\over(1-z_{2})^{2}}+{(\theta_{1}^{2}-\theta_{2}^{2})^{2}\over\langle
w,z\rangle^{2}}\Big{)}>\Big{(}{(\theta_{1}^{1}-\theta_{2}^{1})(\theta_{1}^{2}-\theta_{2}^{2})\over\langle
w,z\rangle^{2}}-{w^{*}\over\langle w,z\rangle}\Big{)}^{2}.$

Observe: $\langle w,z\rangle\geq\min_{i,j}w(i,j)>0$ and so

$\big{|}{w^{*}\over\langle
w,z\rangle}\big{|}\leq{|w^{*}|\over\min_{i,j}w(i,j)}.$

On the other hand, for any $(z_{1},z_{2})\in[0,1]\times[0,1]$,

${\alpha_{1}\over
z_{1}^{2}}+{\alpha_{2}\over(1-z_{1})^{2}}\geq(\alpha_{1}^{1\over
3}+\alpha_{2}^{1\over 3})^{3},\quad{\beta_{1}\over
z_{2}^{2}}+{\beta_{2}\over(1-z_{2})^{2}}\geq(\beta_{1}^{1\over
3}+\beta_{2}^{1\over 3})^{3}.$

Therefore, it is not hard to see that when

$(\alpha_{1}^{1\over 3}+\alpha_{2}^{1\over 3})^{3}(\beta_{1}^{1\over
3}+\beta_{2}^{1\over 3})^{3}>\big{(}{w^{*}\over\min_{i,j}w(i,j)}\big{)}^{2},$
(6.11)

then the Hessian is always positive definite, i.e. the function $g$ in (6.10)
is strictly convex, and the minimum unique. In particular, the condition holds
if $w^{*}=0$. In particular, if $\alpha_{1}=\beta_{1}$, $\alpha_{2}=\beta_{2}$
and $w(1,2)=w(2,1)$, then under (6.11) the unique solution $z_{1},z_{2}$
satisfies $z_{1}=z_{2}$ (by symmetry). Indeed, in symmetric case it holds: if
$(z_{1},z_{2})$ is a solution, then so must be $(z_{2},z_{1})$, and if the
solution is unique, then it must be that $z_{1}=z_{2}$. If (6.11) fails, then
the function might be non-convex, and for small $\alpha_{k}$ and $\beta_{k}$
values it is (given $w^{*}\neq 0$), but evaluated at the minimums, the Hessian
might still be positive definite and so Theorem 6.2 might apply.

Similarly, for (6.8)

$g(z)=-\ln f(z)=\langle\phi,z\rangle-\alpha_{1}^{1}\ln
z_{1}-\alpha_{2}^{1}\ln(1-z_{1})-\alpha_{1}^{2}\ln
z_{2}-\alpha_{2}^{2}\ln(1-z_{2})$ (6.12)

and so the Hessian is

$\displaystyle\left(\begin{array}[]{cc}{\partial^{2}g\over\partial
z_{1}^{2}}&{\partial^{2}g\over\partial z_{1}\partial z_{2}}\\\
{\partial^{2}g\over\partial z_{2}\partial z_{1}}&{\partial^{2}g\over\partial
z_{2}^{2}}\\\ \end{array}\right)=\left(\begin{array}[]{cc}{\alpha_{1}\over
z_{1}^{2}}+{\alpha_{2}\over(1-z_{1})^{2}}&\phi^{*}\\\ \phi^{*}&{\beta\over
z_{2}^{2}}+{\beta_{2}\over(1-z_{2})^{2}}\\\ \end{array}\right),$

where

$\phi^{*}:=\phi(1,1)-\phi(1,2)-\phi(2,1)+\phi(2,2).$

Since the elements of the main diagonal are are strictly positive, the matrix
is positive definite if and only if the determinant is positive:

$\big{(}{\alpha_{1}\over
z_{1}^{2}}+{\alpha_{2}\over(1-z_{1})^{2}}\big{)}\big{(}{\beta_{1}\over
z_{2}^{2}}+{\beta_{2}\over(1-z_{2})^{2}}\big{)}>(\phi^{*})^{2}.$ (6.13)

Thus, when the following inequality holds

$(\alpha_{1}^{1\over 3}+\alpha_{2}^{1\over 3})^{3}(\beta_{1}^{1\over
3}+\beta_{2}^{1\over 3})^{3}>(\phi^{*})^{2},$ (6.14)

then the Hessian is always positive definite and the function (6.12) strictly
convex implying that the minimum is unique. If the Hessian is not always
unique, but (6.13) holds for minimums, then Theorem 6.2 applies. Again, when
$\alpha=\beta$ and $\phi(1,2)=\phi(2,1)$, then under (6.14) the unique minimum
is such that $z_{1}=z_{2}$. For example, when $\alpha=\beta=(2,2)$ and
$\phi(1,1)=1,\phi(1,2)=\phi(2,1)=2,\phi(2,2)=3$, then $\phi^{*}=0$ and,
therefore, (6.14) holds. It means the minimum is unique, $z_{1}=z_{2}$, and
one can verify that

$z_{1}=z_{2}={\sqrt{17}-3\over 2}\approx 0.561.$

So the unique limit distribution in this case is $q\times q$, where
$q=(z,1-z)$. But when $\alpha=(2,3)$, $\beta=(3,2)$ and $\phi$ is as
previously, then the minimum is again unique but since $\alpha\neq\beta$, we
now have $z_{1}\neq z_{2}$: $z_{1}=\sqrt{6}-2\approx 0.45,\quad
z_{2}=\sqrt{7}-2\approx 0.645.$ This means: the unique limit distribution is
$q^{1}\times q^{2}$, where $q^{1}=(\sqrt{6}-2,3-\sqrt{6})$,
$q^{2}=(\sqrt{7}-2,3-\sqrt{7})$.

When $\alpha=\beta=(2,2)$, $\phi(1,1)=4,\phi(1,2)=\phi(2,1)=2,\phi(2,2)=4$,
then $\phi^{*}=4$, but (6.14) still holds and therefore there is unique
minimum: $z_{1}=0.5,z_{2}=0.5$. However, when the $\phi$ is as previously, but
$\alpha=\beta=(0.25,0.25)$, then (6.14) fails. It turns out that now the
function is not convex and there are two minima:

$(z_{1}={2-\sqrt{2}\over 4},z_{2}={2+\sqrt{2}\over
4}),\quad(z_{1}={\sqrt{2}+2\over 4},z_{2}={\sqrt{2}-2\over 4})$

Observe that $z_{2}=1-z_{1}$. Thus the limit measures are $q^{1}\times q^{2}$
and $q^{2}\times q^{1}$, where $q^{1}=(z_{1},z_{2})$ and
$q^{2}=(z_{2},z_{1})$. These two product measures are different. Finally
observe that in both cases (6.13) holds, so that by Theorem 6.2,

$\nu_{n}\Rightarrow{1\over 2}\delta_{q^{1}\times q^{2}}+{1\over
2}\delta_{q^{2}\times q^{1}}.$

The function

$f_{n}(q)=\langle e^{-{\phi\over
n^{\lambda}}},q\rangle^{n^{\lambda}}\prod_{l=1}^{2}\prod_{k=1}^{2}(q(k)^{l})^{0.25-{1\over
n^{1-\lambda}}}$

is symmetric, i.e. $f_{n}(z_{1},z_{2})=f_{n}(z_{2},z_{1})$ and then
$\bar{\pi}_{n}\Rightarrow{1\over 2}\delta_{q^{1}\times q^{2}}+{1\over
2}\delta_{q^{2}\times q^{1}}$. Since now $r_{q}=q$, by Lemma 4.1, (4.8) holds,
with

$P^{*}(x_{1},\ldots,x_{m})={1\over 2}\prod_{i=1}^{m}q^{1}\times
q^{2}(x_{i})+{1\over 2}\prod_{i=1}^{m}q^{2}\times q^{1}(x_{i}),\quad
x_{i}\in{\cal X}^{L}.$

##  7 Experiments

Let $K=2$, $\phi(1)=0$, $\phi(2)=\ln 6$, $\alpha_{1}=0.3$, $\alpha_{2}=0.7$.
Let us find the limit measures $q^{*}$ as in Theorem 6.1 in the following
cases: $\lambda=0$, $\lambda\in(0,1)$ and $\lambda=1$.

##### Case $\lambda=0$:

Then, as it can be easily checked by verifying (6.4) that $q^{*}=({3/5},2/5)$.
Since $\theta=\langle q^{*},w\rangle=2/3$, the measure $r^{*}$ is as follows:
$r^{*}_{1}=w(1)q^{*}(1)/\theta=9/10$ and $r^{*}(2)=w(2)q^{*}(2)/\theta=1/10$.
Therefore the limit process governed by $P^{*}$ is an i.i.d. process with
measure $r^{*}$, and so due to the weight function, the proportion of the
first genotype has increased from 0.3 (according to the prior) to 0.9. Figure
1 illustrates the convergence.

Figure 1: Density histogram of the fraction of the first type in the
population, for a Dirichlet prior $\alpha=(0.3,0.7)$, fitness $\phi=(0,\ln
6)$, and $\lambda=0$, and three population sizes $10^{2}$, $10^{3}$, and
$10^{4}$. The histograms were constructed by recording the fraction of the
first type in the population over $10^{8}$ MCMC samples according to the
process described in section 3.2

.

##### Case $\lambda\in(0,1)$:

The solution of (6.5) is

$q^{*}(1)={\ln(6)-1+\sqrt{(\ln(6)-1)^{2}+1.2\ln(6)}\over 2\ln(6)}\approx
0.686.$ (7.1)

Now $r^{*}=q^{*}$, thus we see that the limit process governed by $P^{*}$ is
an i.i.d. process with measure $q^{*}$, and so due to the weight function, the
proportion of first genotype has increased from 0.3 (according to the prior)
to 0.689. The increase is smaller than in the previous case. Figures 4, 4 and
4 illustrates the convergence for $\lambda=0.25,0.5,0.75$, respectively. We
see that although the limit is the same, the speed of convergence depends very
much on $\lambda$.

Figure 2: Case $\lambda=0.25$ : density histogram of the fraction of the first
type in the population, for a Dirichlet prior $\alpha=(0.3,0.7)$, fitness
$\phi=(0,\ln 6)$, and $\lambda=0.25$, and population sizes $10^{3}$ and
$10^{4}$. The histograms were constructed by recording the fraction of the
fitter allele in the population over $10^{8}$ MCMC samples according to the
process described in section 3.2

.

Figure 3: Case $\lambda=0.5$ : density histogram of the fraction of the first
type in the population, for a Dirichlet prior $\alpha=(0.3,0.7)$, fitness
$\phi=(0,\ln 6)$, and population sizes $10^{3}$ and $10^{4}$. The histograms
were constructed by recording the fraction of the fitter allele in the
population over $10^{8}$ MCMC samples according to the process described in
section 3.2

.

Figure 4: Case $\lambda=0.75$ : density histogram of the fraction of the first
type in the population, for a Dirichlet prior $\alpha=(0.3,0.7)$, fitness
$\phi=(0,\ln 6)$, and population sizes $10^{3}$ and $10^{4}$. The histograms
were constructed by recording the fraction of the fitter allele in the
population over $10^{10}$ MCMC samples according to the process described in
section 3.2

.

Let $q^{*}_{m}=(q_{m}(1),q_{m}(2))$ be the maximizer of

$\ln\langle\exp[-{\phi\over m}],q\big{\rangle}+\sum_{k}{\alpha_{k}\over m}\ln
q(k).$

In our example $q_{1}(1)=3/5$ (the case $\lambda=0$) and

$q_{m}(1)={b_{m}+\sqrt{b_{m}^{2}+{12\over 10}(1+m)\big{(}1-(1/6)^{1\over
m}\big{)}(1/6)^{1\over m}}\over 2(1+m)\big{(}1-(1/6)^{1\over m}\big{)}},\text{
where }b_{m}=(m+0.3)-(1/6)^{1\over m}(1.3+m).$

According to Proposition 6.1, in the process $m\to\infty$, $q_{m}(1)$ tends to
$q^{*}(1)$ from (7.1), and one can easily verify that it is really so.

##### Case $\lambda=1$:

According to 2) of Theorem 5.1, the limit measure $\bar{\pi}$ has density with
respect to the Lebesgue’i measure on $[0,1]$:

$\bar{\pi}(q)={1\over Z}\exp[-\ln 6\cdot q](1-q)^{0.3-1}q^{0.7-1},$

where $Z$ is the valuer of the moment generating function of
Beta(0.7,0.3)-distributed random variable evaluated at $\ln 6$. In the density
above, $q$ stands for $q(2)$. Therefore, the limit proportion of the second
genotype in the stochastic process governed by $P^{*}$ is a random variable,
its distribution has density $\bar{\pi}(q)$ as stated above. Figure 5
illustrates the convergence.

Figure 5: Case $\lambda=1$ : Density histogram (log scale) of the fraction of
the first type in the population, for a Dirichlet prior $\alpha=(0.3,0.7)$,
fitness $\phi=(0,\ln 6)$, and population sizes $10^{2}$ and $10^{3}$. Note
that because the concentration parameter is kept constant, the limit
distribution is a reweighted $\beta$ distribution, with infinities at 0 and 1:
when $\lambda=0$, the prior is the same for all $N$, and the log-fitness term
is reduced by a factor of $N$, so that $P_{n}$ never concentrates. The
histograms were constructed by recording the fraction of the fitter allele in
the population over $10^{9}$ and $10^{10}$ MCMC samples respectively according
to the process described in section 3.2

.

##  8 Conclusions

We have constructed MCMC algorithms that are similar to existing genetic
algorithms. The ‘breeding’ consists of sampling from exchangeable
distributions based on the Dirichlet distribution, and the ‘selection’ is
essentially Metropolis-Hastings. The sequence of populations forms a
reversible Markov chain that satisfies detailed balance conditions. We have
exhibited two possible sampling distributions: more elaborate exchangeable
sampling distributions are possible. The entire MCMC procedure is a population
generalisation of Metropolis-Hastings. As far as we are aware, this is the
first plausible and general model of sexual reproduction that exactly
satisfies detailed balance, and for which the stationary distribution can be
written in closed form for arbitrary fitness functions.

We also explored some properties of the stationary distribution, and showed
that for any fitness function there are three non-trivial limiting
distributions for large population sizes, with two phase transitions. This is
a first step towards a more general understanding of the interaction of the
population size, fitness scaling, and mutation rate in genetic algorithms and
evolutionary models.

Formulating a genetic model as a MCMC procedure opens a new research direction
in using the many techniques developed in MCMC to achieve faster convergence
to the stationary distribution using different MCMC kernels.

We have shown that the stationary distribution is unaffected by multiplicative
noise in fitness evaluations. This has been suggested by, for example, [13],
but our techniques allow a proof of this effect.

Finally there is a more general conclusion from our analysis. For many years,
since [11] and [8], a widely suggested folk-motivation for genetic algorithms
has been that because they are inspired by natural biological evolution, and
because evolution has produced the variety of life on earth, genetic
algorithms should be in some sense generally effective. Our analysis makes it
clear that genetic algorithms are more closely related to conventional MCMC
methods for non-parametric Bayesian inference than has previously been
recognised.

## Appendix

##### Proof of claim 4.1

Proof. First note that

$\sup_{q\in{\cal P}}|\langle w_{n},q\rangle-\langle
w,q\rangle|\leq\|w_{n}-w\|\|q\|\leq\|w_{n}-w\|\to 0.$

Now use the fact that if $f_{n},f,g_{n},g$ are nonnegative functions such that
$\sup_{x}|f_{n}(x)-f(x)|\to 0$, $\sup_{x}|g_{n}(x)-g(x)|\to 0$,
$\inf_{x}g(x)=g_{*}>0$ and $f$, $g$ are bounded above, then with
$g^{*}=\sup_{x}g(x)$ and $f^{*}=\sup_{x}f(x)$

$\displaystyle\sup_{x}\big{|}{f_{n}(x)\over f(x)}-{g_{n}(x)\over g(x)}\big{|}$
$\displaystyle=\sup_{x}\big{|}{f_{n}(x)g(x)-g_{n}(x)f(x)\over
g_{n}(x)g(x)}\big{|}$
$\displaystyle\leq\sup_{x}\big{|}{(f_{n}(x)-f(x))g(x)\over
g_{n}(x)g(x)}\big{|}+\sup_{x}\big{|}{(g_{n}(x)-g(x))f(x)\over
g_{n}(x)g(x)}\big{|}$
$\displaystyle\leq\sup_{x}\big{|}{(f_{n}(x)-f(x))g^{*}\over
g_{n}(x)g(x)}\big{|}+\sup_{x}\big{|}{(g_{n}(x)-g(x))f^{*}\over
g_{n}(x)g(x)}\big{|}\to 0,$

because for $n$ big enough $g_{n}(x)g(x)\geq{g^{2}_{*}\over 2}$ for every $x$.
Take $x=q$, $f_{n}(q)=w_{n}(k)q(k)$, $f(q)=w(k)q(k)$, $g_{n}(q)=\langle
w_{n},q\rangle$ and $g(q)=\langle w,q\rangle$. Then $g_{*}=w(K)>0$,
$f^{*}=g^{*}=w(1)$ and so (4.9) follows.

###  8.1 Proof of Proposition 5.1

Recall that $f_{n}$ and $f$ are continuous and bounded functions on ${\cal P}$
so that $\|f_{n}\|_{\infty}<\infty$ and $\|f\|_{\infty}<\infty$. By
assumption, $\pi$ is a finite measure. Since $f_{n}$ converges to $f$
uniformly, it follows that $\|f_{n}-f\|_{\infty}\to 0$ and so
$\|f_{n}\|_{\infty}\to\|f\|_{\infty}.$ For every $m$,

$|\|f_{n}\|_{m}-\|f\|_{m}|\leq\|f_{n}-f\|_{m}\leq\pi({\cal P})^{1\over
m}\|f_{n}-f\|_{\infty}\to 0.$

Since $\|f\|_{m_{n}}\to\|f\|_{\infty}$, we have

$\displaystyle\big{|}\|f_{n}\|_{m_{n}}-\|f\|_{\infty}\big{|}$
$\displaystyle\leq\big{|}\|f_{n}\|_{m_{n}}-\|f\|_{m_{n}}\big{|}+\big{|}\|f\|_{m_{n}}-|f\|_{\infty}\big{|}$
$\displaystyle\leq\|f_{n}-f\|_{m_{n}}+\big{|}\|f\|_{m_{n}}-|f\|_{\infty}\big{|}\leq\pi({\cal
P})^{1\over
m_{n}}\|f_{n}-f\|_{\infty}+\big{|}\|f\|_{m_{n}}-|f\|_{\infty}\big{|}\to 0.$

Now fix $\delta>0$. Since ${\cal
S}^{*}_{\delta}:=\\{q:f(q)>\|f\|_{\infty}-\delta\\}$, we have ${\cal S}-{\cal
S}^{*}_{\delta}=\\{q:f(q)\leq\|f\|_{\infty}-\delta\\}.$ Define
$\delta^{\prime}:=\delta/\|f\|_{\infty}$. Then

$\displaystyle\sup_{q\in{\cal S}-{\cal
S}^{*}_{\delta}}{f_{n}(q)\over\|f_{n}\|_{m_{n}}}$
$\displaystyle=\sup_{q\in{\cal S}-{\cal
S}^{*}_{\delta}}{f(q)+(f_{n}(q)-f(q))\over\|f_{n}\|_{m_{n}}}=\sup_{q\in{\cal
S}-{\cal
S}^{*}_{\delta}}{f(q)+(f_{n}(q)-f(q))\over\|f\|_{\infty}}{\|f\|_{\infty}\over\|f_{n}\|_{m_{n}}}$
$\displaystyle\leq\sup_{q\in{\cal S}-{\cal
S}^{*}_{\delta}}{f(q)\over\|f\|_{\infty}}{\|f\|_{\infty}\over\|f_{n}\|_{m_{n}}}+{\|f_{n}-f\|_{\infty}\over\|f_{n}\|_{m_{n}}}\leq
1-{\delta^{\prime}\over 2},$

provided $n$ is big enough. Thus,

$\sup_{q\in{\cal S}-{\cal
S}^{*}_{\delta}}h_{n}(q)\leq\Big{(}1-{\delta^{\prime}\over
2}\Big{)}^{m_{n}}\to 0,$

so that $\nu_{n}({\cal S}^{*}_{\delta})\to 1$. We now argue that for any
$\epsilon>0$ there exists $\delta>0$ so that

${\cal S}^{*}_{\delta}\subset B(q^{*},\epsilon),$ (8.1)

where $B(q^{*},\epsilon)$ is a ball in Euclidean sense. If, for an
$\epsilon>0$, such a $\delta>0$ would not exists, then there would exist a
sequence $q_{n}\to q$ such that $f(q_{n})\nearrow f(q)$, but
$\|q_{n}-q\|\geq\epsilon.$ Since ${\cal P}$ is compact, along a subsequence
$q_{n^{\prime}}\to q$ and by continuity $f(q_{n^{\prime}})\to f(q)$. On the
other hand $\|q-q^{*}\|>\epsilon$ and that would contradict the uniqueness of
$q^{*}$. Therefore (8.1) holds and so for any $\epsilon>0$, it holds that
$\nu_{n}\big{(}B(q^{*},\epsilon)\big{)}\to 1$. From the definition of the weak
convergence, it now follows that $\nu_{n}\Rightarrow\delta_{q^{*}}.$

###  8.2 Proof of Lemma 6.1

1)

To find

$q^{*}=\arg\max_{q\in{\cal P}}\,\,[\ln\langle
e^{-\phi},q\rangle+\sum_{k}\alpha_{k}\ln q(k)],$ (8.2)

we define Lagrangian

$L(q,\beta)=\ln\langle e^{-\phi},q\rangle+\sum_{k}\alpha_{k}\ln
q(k)-\beta(\sum_{k}q(k)-1)$

(here $\beta$ is a scalar) and maximize $L(q,\beta)$ over $q>0$ (all entries
are positive). Taking partial derivatives with respect to $q(k)$, we have

${e^{-\phi(k)}\over\langle e^{-\phi},q\rangle}+{\alpha_{k}\over
q(k)}=\beta,\quad\Rightarrow\quad{e^{-\phi(k)}q(k)\over\langle
e^{-\phi},q\rangle}+{\alpha_{k}}=q(k)\beta\quad\forall k.$

With $|\alpha|=\sum_{k}\alpha_{k},$ we have thus $\beta=1+|\alpha|$ and so the
solution $q^{*}$ satisfies the set of equalities

$q^{*}(k)={1\over 1+|\alpha|}\Big{(}{e^{-\phi(k)}q^{*}_{k}\over\langle
e^{-\phi},q^{*}\rangle}+\alpha_{k}\Big{)},\quad\forall k.$ (8.3)

Now with $w(k)=e^{-\phi(k)}$ define parameter $\theta:=\langle w,q^{*}\rangle$
and rewrite (8.3) as follows

$q^{*}(k)={\alpha_{k}\over(1+|\alpha|)-{w(k)\over\theta}}\quad k=1,\ldots,K.$
(8.4)

We see that amongst the probability vectors satisfying $\langle
w,q^{*}\rangle=\theta,$ the solution is unique. Since $\alpha_{k}>0$ for every
$k$, it is easy to see that there is only one parameter $\theta$ such that the
right hand side of (8.4) would be a probability measure: if
$\theta^{\prime}>\theta$, then for every $k$, we have

${\alpha_{k}\over(1+|\alpha|)-{w(k)\over\theta}}>{\alpha_{k}\over(1+|\alpha|)-{w(k)\over\theta^{\prime}}}.$

Therefore a solution of (8.2) is unique vector $q^{*}$ given by (8.4), where
$\theta=\langle w,q^{*}\rangle$.

2)

To find

$q^{*}=\arg\max_{q\in{\cal P}}[-\langle\phi,q\rangle+\sum_{k}\alpha_{k}\ln
q(k)],$ (8.5)

we define Lagrangian

$L(q,\beta)=-\langle\phi,q\rangle+\sum_{k}\alpha_{k}\ln
q(k)-\beta(\sum_{k}q(k)-1).$

Partial derivatives with respect to $q(k)$ give us the equalities

$-\phi(k)+{\alpha_{k}\over q(k)}=\beta\quad\forall
k\quad\Rightarrow\quad-\langle\phi,q\rangle+|\alpha|=\beta.$

Therefore, the inequalities for $q^{*}(k)$ are

$q^{*}(k)={\phi(k)q^{*}(k)-\alpha_{k}\over\langle\phi,q^{*}\rangle-|\alpha|}={\phi(k)q^{*}(k)-\alpha_{k}\over\theta-|\alpha|},\quad\theta:=\langle\phi,q^{*}\rangle.$
(8.6)

After rewriting (8.6), we obtain

$q^{*}(k)={\alpha_{k}\over\phi(k)+|\alpha|-\theta},\quad k=1,\ldots,K,$

Thus, there cannot be two solutions having the same $\theta$. As in the case
1), it is easy to see that when $\alpha_{k}>0$ there is only one $\theta$ so
that (2)) sums up to one. Therefore, the solution to the problem (8.5) is
unique. Note that the solution is independent of $\lambda$.

### Acknowledgments

The research is supported by Estonian Institutional research funding IUT34-5.

## References

  * [1] S. Baluja and R. Caruana.  Removing the genetics from the standard genetic algorithm.  In ICML, pages 38–46. Morgan Kauffman Publishers, Inc., 1995.
  * [2] Shumeet Baluja.  Genetic algorithms and explicit search statistics.  In Advances in Neural Information Processing Systems, pages 319–325, 1997.
  * [3] E.B. Baum, D. Boneh, and C. Garrett.  Where genetic algorithms excel.  Evolutionary Computation, 9(1):93–124, 2001.
  * [4] James Franklin Crow and Motoo Kimura.  An introduction to population genetics theory. New York, Evanston and London: Harper & Row, Publishers, 1970.
  * [5] Lloyd Elliott and Yee Whye Teh.  Scalable imputation of genetic data with a discrete fragmentation-coagulation process.  In Advances in Neural Information Processing Systems 25, pages 2861–2869, 2012.
  * [6] W.J. Ewens.  Mathematical population genetics: theoretical introduction, volume 1.  Springer Verlag, 2004.
  * [7] Stuart Geman and Donald Geman.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  Pattern Analysis and Machine Intelligence, IEEE Transactions on, (6):721–741, 1984.
  * [8] D.E. Goldberg.  Genetic algorithms in search, optimization, and machine learning.  Addison-wesley, 1989.
  * [9] Heikki Haario and Eero Saksman.  Simulated annealing process in general state space.  Advances in Applied Probability, 23(4):866–893, 1991.
  * [10] Nils Lid Hjort, Chris Holmes, Peter Muller, and Stephen Walker.  Bayesian nonparametrics.  Number 28. Cambridge University Press, 2010.
  * [11] J.H. Holland.  Adaptation in natural and artificial systems.  University of Michigan press, 1975.
  * [12] P.A.P. Moran.  The statistical processes of evolutionary theory. Clarendon Press; Oxford University Press., 1962.
  * [13] Gregory Morse and Kenneth O Stanley.  Simple evolutionary optimization can rival stochastic gradient descent in neural networks.  In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pages 477–484. ACM, 2016.
  * [14] R.M. Neal.  Markov chain sampling methods for Dirichlet process mixture models.  Journal of computational and graphical statistics, pages 249–265, 2000.
  * [15] M.J.A. Strens.  Evolutionary MCMC sampling and optimization in discrete spaces.  In ICML, volume 20, page 736, 2003.
  * [16] Yee Whye Teh and Michael I Jordan.  Hierarchical Bayesian nonparametric models with applications.  In Bayesian nonparametrics, pages 158–207. Cambridge University Press, 2010.
  * [17] Cajo JF Ter Braak.  A Markov Chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces.  Statistics and Computing, 16(3):239–249, 2006.
  * [18] M.D. Vose.  The simple genetic algorithm: foundations and theory.  The MIT Press, 1999.
  * [19] Chris Watkins and Yvonne Buttkewitz.  Sex as Gibbs sampling: a probability model of evolution.  arXiv preprint arXiv:1402.2704, 2014.