An evolutionary model that satisfies detailed balance

TL;DR

This paper introduces a class of evolutionary models with exchangeable breeding and fitness-based selection, enabling explicit analysis of mutation-selection equilibrium and phase transitions as population size grows.

Contribution

It presents a novel evolutionary model framework where the stationary distribution can be explicitly derived and studied, including phase transition phenomena.

Findings

Explicit stationary distribution derived

Phase transition theorems proved

Behavior analyzed as population size increases

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.

Equations

$$P_{\xi }(\mathbf{x}):=P_{\xi }(x_1,\dots ,x_n):=P({\xi }_1=x_1,\dots ,{\xi }_n=x_n)$$

$$P_{\xi }(x_1,\dots ,x_n)=P_{\xi }(x_{\sigma (1)},\dots ,x_{\sigma (n)}).$$

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

$$\underset{\text{stationary distribution}}{\underbrace{P_n(x_1,\dots ,x_n)}}=\frac{1}{Z_n}\underset{\text{breeding term}}{\underbrace{P_{\xi }(x_1,\dots ,x_n)}}\underset{\text{fitness term}}{\underbrace{w(x_1)\cdots w(x_n)}}.$$

$$\mathcal{P}:=\{(q_1,\dots ,q_k):q_i\ge 0,\sum _i{q}_i=1\},$$

$$P_{\xi }(\mathbf{x})={\int }_{\mathcal{P}}\prod _{i=1}^{n}q(x_i)\pi (dq).$$

$$(X_{1,n},\dots X_{m,n})\Rightarrow (X_1,\dots X_m)?$$

$$w_n(k)=\exp \left[-\frac{\varphi (k)}{n^{\lambda }}\right],k=1,\dots ,K$$

$$\pi (q)\exp [-⟨\varphi ,q⟩].$$

$${\pi }_n=\mathrm{Dir}(n^{1-\lambda }{\alpha }_1,\dots n^{1-\lambda }{\alpha }_K),$$

$$\pi (q_1,\dots ,q_k)\propto \prod _{k=1}^K{q}_k^{n^{1-\lambda }{\alpha }_k}.$$

$$\exp [-⟨\varphi ,q⟩]\prod _{k=1}^K{q}_k^{{\alpha }_k-1}.$$

$$n_k=\sum _{i=1}^n{I}_{\{k\}}(x_i).$$

$$P_{\xi }({\xi }_{n+1}=k\mid x_1,\dots ,x_n,\alpha ):=\frac{n_k+{\alpha }_k}{n+|\alpha |}.$$

$$P_{\xi }(x_1,\dots ,x_n\mid \alpha )=\frac{({\alpha }_1{)}_{n_1}\cdots ({\alpha }_K{)}_{n_K}}{(|\alpha |{)}_n},$$

$$(\alpha {)}_n:=\alpha (\alpha +1)\cdots (\alpha +n-1).$$

$$P_{\xi }(n_1,\dots ,n_K)=\frac{n!}{n_1!\cdots n_K!}\frac{({\alpha }_1{)}_{n_1}\cdots ({\alpha }_K{)}_{n_K}}{(|\alpha |{)}_n}.$$

$$P_{\xi }(x_1,\dots ,x_n\mid \alpha ):={\int }_{\mathcal{P}}q(x_1)\cdots q(x_n)\pi (dq),$$

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

$$u_k:=\frac{{\alpha }_k}{|\alpha |+n}\text{or equivalently}{\alpha }_i=\frac{n{u}_i}{1-u}.$$

$$\frac{n_i+{\alpha }_i}{|\alpha |+n}=\frac{n_i}{n}\frac{n}{|\alpha |+n}+\frac{{\alpha }_i}{|\alpha |+n_i}=\frac{n_i}{n}(1-u)+u_i$$

$$\frac{n_i{n}_j}{n^2}(1-u)+\frac{n_j}{n}{u}_i.$$

$$\frac{n_1{n}_2}{n^2}(1-u_2)+\frac{n_2^2}{n^2}{u}_1,$$

$$\frac{n_i{n}_j}{n^2}+u_i\frac{n_j}{n}$$

$$P_{\xi }(x_1,\dots ,x_n):=\prod _{j=1}^L{P}_{{\xi }^j}(y_1^j,\dots ,y_n^j)$$

$$\mathbf{x}$$

$$P_{\xi }(\mathbf{x})P_{\xi }(x_{n+1}\mid \mathbf{x})=P_{\xi }(x_1,\dots ,x_{n+1})=P_{\xi }({\mathbf{x}}^{\prime })P_{\xi }(x_i\mid {\mathbf{x}}^{\prime })$$

$$P_n(\mathbf{x})P(\mathbf{x}\to {\mathbf{x}}^{\prime })$$

$$=\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}}$$

$$=\frac{1}{Z_n}{P}_{\xi }(x_1,\dots ,x_{n+1})\cdot \frac{1}{n}\frac{w_1\cdots w_{n+1}}{w_i+w_{n+1}}$$

$$=\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}}$$

$$=\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}}$$

Full text

Keywords:

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

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 analyze 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 stationary distribution is.

For Moran models with only a single locus, the stationary distribution can be characterised. But for models with more than one locus, and fitness that depends nonlinearly on the combination of alleles at the loci, the stationary distribution is notoriously hard to characterise. The reason is that in previous models of sexual reproduction with multiple loci, 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 [25], but these calculations are neither easy nor revealing for arbitrary fitness functions on multiple loci.

For a reversible Markov chain, the stationary distribution may be found by verifying detailed balance conditions. In our model, introduced in [26], 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 kind of evolution is typically modeled with a Moran model [18] and in subsections 2.1 and 3.4, we shall see that when the breeding process is a generalised Polya urn model, then our model can be indeed considered as a discrete time Moran model with mutation and selection. But our model is more general that just a Moran one, because we allow a more general breeding process. However, due to the connection with Moran models, the generalized Polya urn model (or, equivalently, Dirichlet-categorical process) is an important special case of breeding so that section 6 is devoted to that model.

In many evolutionary models such as [5], and in genetic algorithms such as [15], the reproductive fitness of a genome is modelled as the genome’s rate of breeding, and not as its probability of death. In these models, in the breeding phase of each generation, fit genomes are chosen more frequently to breed than unfit genomes are: in these models, discarding individuals – or ‘death’ – is modelled as uniformly random deletion from the population. The second option is that the selection occurs at death: individuals are discarded (or ‘die’) with probability related to their fitness (termed viability selection): less-fit individuals are more likely to be discarded. For discrete time classical models such as the Wright-Fisher or Moran model, there is actually no big difference whether the selections occur at death or at reproduction (although formulas can be slightly different, see e.g. [20]) and in the literature both versions are present. In our model, breeding is conditional sampling from a very general exchangeable process and so it is quite natural that selection occurs at death. This is a significant design choice that turns out to greatly simplify the analysis. Observe that selection at death it does not mean that all allele types have equal probability to be born. Our model might have breeding process as i.i.d. random variables with given distribution possibly very far form uniform, or the breeding process can be generalized Polya urn, where the history of evolution determines what types are more likely to be born.

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 Metropolis-Hastings process.

Mathematical formulation, overview of the main results and

organization of the paper.

Throughout the paper, we suppose there is a finite set ${\cal X}=\{1,\ldots,K\}$ of possible genomes. Our models apply both to the case of genomes with one locus with $K$ possible alleles, and also to the case where genomes have more than one locus: for example, for genomes with $L$ loci with two possible alleles at each, then $K=2^{L}$. Throughout the paper, $\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

[TABLE]

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

[TABLE]

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:

[TABLE]

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. As explained in subsections 3.1, 3.2 and 3.3, there are many possible schemes to formalize this idea, and we shall show that all those schemes lead to one particular stationary distribution of populations that factorises into the form:

[TABLE]

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

The measure $P_{n}$ is the main object of interest. Since $\xi$ is exchangeable, there exists a prior measure $\pi$ on the set of all probability measures on ${\cal X}$ (called simplex), denoted by

[TABLE]

such that for every $\mathbf{x}\in{\cal X}^{n}$

[TABLE]

Here $q(k)=q_{k}$, in what follows, we shall use both notations ($q_{k}$ and $q(k)$). Hence $P_{n}$ is fully determined by prior measure $\pi$ and fitness function $w$. In order to achieve more generality, we allow $w$ and $\pi$ to depend on population size $n$ (thus writing $w_{n}$ and $\pi_{n}$), and we ask: do the measures $P_{n}$ converge (in some sense) to a limit? The sense of convergence needs to be specified, because $P_{n}$ is defined on ${\cal X}^{n}$, so that the domain of the measure depends on $n$. Therefore we consider the array of random variables $(X_{n,1},\ldots,X_{n,n})\sim P_{n}$ and we ask: is there a stochastic process $X_{1},X_{2},\ldots$ so that for every $m$

[TABLE]

Another option to define the limit is based on the observation that the measure $P_{n}$ is invariant with respect to the permutations (finitely exchangeable), and so it depends on the counts of vector ${\bf x}$, only. This, in turn, allows to define a counterpart of $P_{n}$ – a probability measure $Q_{n}$ on the simplex ${\cal P}$. The formal definition of $Q_{n}$ is given in subsection 4.1. Since all measures $Q_{n}$ are defined on the same domain (simplex ${\cal P}$), we can now ask whether the sequence $Q_{n}$ converges in the classical sense of weak convergence of probability measures. Throughout the paper, these two approaches – measures $P_{n}$ (random process) and measures $Q_{n}$ – are handled in parallel.

Our first limit result, Theorem 5.1 considers the case when the prior measure $\pi$ is independent of $n$, but arbitrary; the fitness function $w_{n}$ depends on $n$ in the following way (recall ${\cal X}=\{1,\ldots,K\}$):

[TABLE]

where $\lambda\geq 0$ and $0\leq\phi(1)\leq\phi(2)\leq\cdots\leq\phi(K)$. The case $\lambda=0$ corresponds to fitness that is constant in that it does not vary with $n$. Observe that for any $n$, the genotype 1 is the fittest. Theorem 5.1 shows that the following phase transition occurs with respect to $\lambda$:

when $\lambda\in[0,1)$ and $\phi(1)<\phi(2)$, then the limit process $X_{1},X_{2},\ldots$ has one trajectory $1,1,1,\ldots$ a.s., and $Q_{n}\Rightarrow\delta_{q^{*}}$, where $q^{*}=(1,\ldots,0).$ Thus, when $\lambda\in[0,1)$ only the fittest genotype survives in the limit, no matter what the prior says; the influence of the prior vanishes and fitness prevails. The additional assumption $\phi(1)<\phi(2)$ means that $w_{n}(1)>w_{n}(k)$, $k=2,\ldots,K$, and so it is quite natural for the result. 2. 2.

When $\lambda=1$, $Q_{n}$ converges to a nondegenerate distribution specified in Theorem 5.1 that depends on the prior $\pi$ and on the function $\phi$. When the prior $\pi$ has density $\pi(q)$, then the limit distribution also has density proportional to

[TABLE]

Then also the limit process $X_{1},X_{2},\ldots$ is exchangeable process with prior (1.5). non-degenerate. 3. 3.

When $\lambda>1$, then the influence of fitness vanishes and only the prior matters – the limit process is equal to the breeding process and $Q_{n}\Rightarrow\pi$.

Section 6 considers Dirichlet priors

[TABLE]

where $\alpha:=(\alpha_{1},\ldots,\alpha_{K})$, $\alpha_{k}>0$. Thus $\pi_{n}$ admits a density on ${\cal P}$

[TABLE]

With a Dirichlet prior the breeding process is a (generalized) Polya urn process, also known as a Dirichlet-categorical process. This process is a natural choice because it admits a notion of ‘mutation’. Conditional sampling from it is directly interpretable as a simplified model of sexual breeding with mutation, as explained in section 2. The fitness is previously defined as in (1.4) and the common parameter $\lambda$ allows to interpolate between influence of fitness and prior. The case $\lambda=1$ corresponds to the constant prior case. Our second main limit theorem, Theorem 6.1 shows that with respect to $\lambda$, phase transition occurs again:

when $\lambda=0$, then the limit process $X_{1},X_{2},\ldots$ consists of i.i.d. random variables with certain distribution $r^{*}$ (specified in Theorem 6.1) depending on $\phi$ and $\alpha$, and $Q_{n}\Rightarrow\delta_{r^{*}}$. 2. 2.

when $\lambda\in(0,1)$, then then the limit process $X_{1},X_{2},\ldots$ consists of i.i.d. random variables with another distribution $q^{*}$ (specified in Theorem 6.1) also depending on $\phi$ and $\alpha$, and $Q_{n}\Rightarrow\delta_{q^{*}}$; 3. 3.

when $\lambda=1$ and $\phi(1)<\phi(2)$, then the limit is specified by Theorem 5.1. From (1.5) it follows that the limit distribution of allele proportions has density proportional to

[TABLE]

The limits of type (1.7) are common for standard models in evolution theory, both for Wright-Fisher and Moran model (see, e.g. [6, 9, 10]). We shall discuss the connection with our and Moran model in subsections 2.1 and 3.4.

The phase transitions at $\lambda=0$ and $\lambda=1$ are sharp. These are the most important results of the paper.

The paper is organized 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 (1.2). Subsection 3.5 establishes that the stationary distributions are invariant to multiplicative fitness noise; subsection 3.6 shows that special cases of this process are forms of Bayesian inference by MCMC. Next, section 4 establishes basic conditions on the convergence of the stationary distribution to a limit distribution as population size tends to infinity. Section 5 is devoted to Theorem 5.1 and in section 6 generalized Polya urn processes (Dirichlet priors) are considered. In this section, Theorem 6.1 is proved and the product of independent Polya urn processes are considered (subsection 6.2). 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 [11] from an exchangeable distribution; exchangeable Gibbs sampling is a standard technique in statistical nonparametric MCMC methods, for example [21, 14], 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 (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.

2.1 Dirichlet-Categorical Process

In the simplest case, each ‘genome’ consists of only one locus or ‘gene’ which can be one of $K$ possible allele ‘types’, 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})$:

[TABLE]

From the Polya urn, the next allele is drawn of type $k$ with the following probability:

[TABLE]

It follows that

[TABLE]

where for any $\alpha>0$

[TABLE]

The process $\xi_{1},\xi_{2},\ldots$ is infinitely exchangeable by inspection. Due to the exchangeability, by the formula above, it is clear that the random vector of counts has the distribution

[TABLE]

By de Finetti’s theorem

[TABLE]

where 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, mutation rate and Moran models.

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 sample (thus $n+|\alpha|$ 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 the $n$ actual balls present: since a new color may be introduced in this way, we regard this as analogous to a mutation. The mutation rate $u$ is

[TABLE]

To be more detailed, for any type $k=1,\ldots,K$, the ‘prior’ ball of type $k$ is sampled with probability

[TABLE]

Hence $u_{1}+\cdots+u_{K}=u$ and $u_{k}$ is the rate (probability) of mutating from any other type to $k$. 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 $u$ but with different values of $n$, $\alpha$ must be adjusted to depend on $n$. Such single birth-death changes are often modeled via the well known Moran model, and it turns out that our model of breeding in some sense generalizes the drift-mutation Moran model. Let us explain that connection more closely. Recall that a standard (simple) $K$-allele Moran model without mutation and selection is as follows. We consider a population of $n$ individuals and suppose that at time $t$, the number of type $k$ individuals is $n_{k}$. An individual is randomly chosen to reproduce and another is independently chosen to die. When the type $i$ is chosen to reproduce and type $j$ is chosen to die, then the allele counts at time $t+1$ are $(n_{1},\ldots,n_{i}+1,\ldots,n_{j}-1,\ldots n_{K})$ and the probability of such a change is ${n_{i}n_{j}\over n^{2}}$. There are several ways to incorporate mutation into the Moran model, but in our model the probability that type $i$ is born is

[TABLE]

In the sum above, the first part is the transition probability without mutation multiplied with the probability of no mutation, the second part is the probability that the selected type mutates to type $i$. When the probability of $j$-type to die is just ${n_{j}\over n}$ we end up with a Moran model, where the probability that $n_{i}$ becomes $n_{i}+1$ and $n_{j}$ becomes $n_{j}-1$ is

[TABLE]

It is easy to verify that a Markov chain with transition (2.5) satisfies the detailed balance equation with stationary distribution being (2.2), where $\alpha_{k}$ is defined as in (2.3). There are surely other ways to specify mutation in Moran model so that (2.2) is the stationary distribution. For $K=2$, one such example is in [9]: let $u_{1}$ be the probability that type $2$ mutates to 1 and $u_{2}$ be the probability that type 1 mutates to 2. Then the transition from $(n_{1},n_{2})$ to $(n_{1}+1,n_{2}-1)$ has probability

[TABLE]

and detailed balance equation with (2.2) (with $\alpha_{k}$ being defined as in (2.3)) is easy to verify (this is (3.58) in [9]). Other examples of Moran models without selection having (2.2) as stationary distribution can be found in [24], where the authors consider a discrete time version of so-called two-allele decoupled Moran model introduced in [4] (see also [22]). In their model, the transition probability from $n_{i}$ and $n_{j}$ to $n_{i}+1$ and $n_{j}-1$ without selection is very close to (2.5):

[TABLE]

and the the detailed balance equation with (2.2) is easy to to verify. In this case, $\alpha_{i}=nu_{i}$. A $K$-allele version of that model with transition matrix as in (2.7) is studied in [8] and the stationary distribution (14) in [8] is exactly (2.2). Hence our breeding model largely incorporates mutation-drift version of Moran models, but it is much general, because the Dirichlet distribution is only one possible choice of prior measure $\pi$ in (1.3).

Note that in our model, mutations only occur at birth, and – importantly – the distribution of mutations does not depend on the frequencies of different colors that are currently in the population; mutations are distributed according to a fixed prior distribution determined by the frequencies of ‘prior’ balls in the urn.

In our model, we count as a mutation any ball which results from a draw of a ‘prior’ 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.2 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

[TABLE]

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 [5, 9].

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 [7]; 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, for brevity 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_{n+1}}{w_{n+1}+w_{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_{n+1}}{w_{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:

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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

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

Sample $i$ from $\{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:

[TABLE]

Note in passing that

[TABLE]

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

[TABLE]

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:

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 The selection schemes and Moran models

Recall that under the Dirichlet prior our breeding process can be considered as Moran model with birth and mutation probabilities given in (2.4). Including the single tournament and inverse fitness selections into the model, we end up with Moran model, where the transition probability that $n_{i}$becomes $n_{i}+1$ and $n_{j}$ becomes $n_{j}-1$ are, respectively,

[TABLE]

By the same arguments as in subsections 3.1 and 3.2, it is easy to verify that both transition matrices satisfy detailed balance equations with the stationary measure being (1.2):

[TABLE]

Also, when adding the same single tournament or inverse fitness selections to the two allele model (2.6), the detailed balance equations with (3.5) hold. The same holds for the decoupled Moran model (2.7): when adding single tournament or inverse fitness selections, the detailed balance equations with (3.5) hold. Of course, our selection schemes are far from the only possibilities for having (3.5) as stationary distribution. For example, when the mutation-drift Moran model (i.e. the probability of $j$-type to eject just being $n_{j}/n$, no selection) satisfies the detailed balance equation with (2.2) being the stationary distribution, then multiplying the transition probability that type $i$ is chosen to reproduce and $j$ to die simply by $w_{i}$ (as it is often done in the literature, e.g. [4, 24, 6]) then the detailed balance equation with (3.5) as stationary distribution holds. Thus in single tournament selection, the term ${w_{i}\over(w_{i}+w_{j})}$ can be replaced by $w_{i}$ and, again, detailed balance equations with (3.5) as stationary distribution hold (independently of drift-mutation probabilities, given that they satisfy detailed balance equation with (2.2)). To summarize: there are various different kinds of Moran models yielding the stationary distribution $P_{n}$ with $\Pi$ being the Dirichlet distribution

Typically in the evolution literature, for haploid models the weights $w_{k}$ are in form $(1+s_{k})$, in two allele models simply $(1+s)$ and 1, where $s$ is either a positive or negative small number. In the large population limit, typically $s\cdot n\to\gamma$, where $\gamma$ is constant. Hence, for large $n$, $(1+s)\approx(1+{\gamma\over n})\approx\exp[{\gamma\over n}]$ and this justifies our choice of weights in (1.4). So, when in two allele model, $\phi(1)=0$ and $\phi(2)=\gamma>0$, thus $w_{n}(1)=1$ and $w_{n}(2)=e^{-\gamma/n}$, then the limit measure (1.7) has density

[TABLE]

that (with negative $\gamma$) is a classical and well-known limit, sometimes called Wright’s formula [4] (see also (9) in [24], Ch 5 in [9], (7.28) in [6], Ch 1 in [10], Ch 4 in [16]). We would like to stress that we obtain our limit without any diffusion approximation and our main theorems 5.1 and 6.1 generalize this result in various ways.

3.5 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 formalize 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}$.

[TABLE]

The joint probability

[TABLE]

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

[TABLE]

Thus the joint measure factorizes

[TABLE]

where

[TABLE]

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

[TABLE]

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.6 Fitness as likelihood: a connection with non-parametric Bayesian MCMC

In Dirichlet Process mixture models, as described for example in [21, 23], 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

[TABLE]

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

[TABLE]

then this distribution (3.6) 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.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}$ in (1.2) 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}$, the measure $P_{\xi}(\mathbf{x})$ allows the representation (1.3). Hence we can write $P_{n}$ as follows

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.4) 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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

We now present a general lemma that guarantees the convergence (4.7). 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

[TABLE]

We start with the following observation, proven in appendix.

Claim 4.1

If $w_{n}(k)\to w(ik)$ $\forall k\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.

[TABLE]

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

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.7) holds. Moreover, for every $(x_{1},\ldots,x_{m})\in{\cal X}^{m}$

[TABLE]

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

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

[TABLE]

Since the functions

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

We may now write:

[TABLE]

In what follows, let us denote

[TABLE]

Observe that

[TABLE]

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

[TABLE]

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)$

[TABLE]

Thus $Q_{n}$ is a discrete measure

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

because

[TABLE]

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.11) and $Q_{n}$ via $P_{n}$ as in (4.12).

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

[TABLE]

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

[TABLE]

where

[TABLE]

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}$:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

By (4.14),

[TABLE]

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

[TABLE]

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$ as defined in (1.4), i.e. $w(k)=\exp(-\phi(k)/n^{\lambda})$. 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.3) with $w_{n}$:

[TABLE]

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}$. When $w(1)>w(2)$, the function $q\mapsto\langle q,w\rangle$ has unique maximizer $q^{*}:=(1,0,\ldots,0)$. The following theorem states that under $w(1)>w(2)$ 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

[TABLE]

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$. When $\lambda=1$, then $\bar{\pi}_{n}$ as well $Q_{n}$ converges to a nondegenerate distribution, when $\lambda>1$, then $\bar{\pi}_{n}\Rightarrow\pi$, $Q_{n}\Rightarrow\pi$ and the measure $P^{*}$ is the law of birth process $\xi$.

Theorem 5.1

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

1)

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

[TABLE]

2)

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

[TABLE]

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

[TABLE]

3)

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

[TABLE]

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}$ having full support (i.e. the support of $\pi$ is ${\cal P}$), and we are interested in the asymptotic behavior of the measure $\nu_{n}$, where

[TABLE]

Here $\|f\|_{m}$ stands for $L_{m}$-norm and 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

[TABLE]

Here $\|f\|_{\infty}$ is the essential supremum of $f$ with respect to the $\pi$-measure. If $f$ is continuous, then $\|f\|_{\infty}=\sup_{q}f(q)$ (because $\pi$ has full support). 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}$ having full support. Then for every $\delta>0$, \nu_{n}\big{(}{\cal P}_{\delta}^{*}\big{)}\to 1. If ${\cal P}^{*}=\{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

[TABLE]

Proof. (Theorem 5.1)

1)

For $\lambda=0$, take $f_{n}(q)=f(q)=\langle w,q\rangle$, $m_{n}=n$. By $w(1)>w(2)\geq\cdots\geq w(K)$, we have ${\cal P}^{*}=\{q^{*}\}$ , 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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Recall

[TABLE]

Therefore, from (5.4), it follows that when $\lambda=1$, we have for every measurable $E$

[TABLE]

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

[TABLE]

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$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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\}$, $D$ stands for Kullback-Leibler divergence and $c>0$ is a constant.

6 Dirichlet prior

Also in the current section we consider the weights $w_{n}(k)$ as in (1.4), 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. Therefore, it is meaningful to consider non-constant priors so that the influence of prior increases with suitable rate. This motivates us to consider Dirichlet priors (1.6), i.e. $\pi_{n}={\rm Dir}(n^{1-\lambda}\alpha_{1},\ldots n^{1-\lambda}\alpha_{K}),$ 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 its 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 (1.4) and the prior $\pi_{n}$ as in (1.6). Let $\bar{\pi}_{n}$ be defined as in (5.1) 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

[TABLE]

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

[TABLE]

2)

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

[TABLE]

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

[TABLE]

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

Lemma 6.1

1)

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

[TABLE]

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

2)

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

[TABLE]

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:

[TABLE]

where

[TABLE]

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

[TABLE]

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.3), i.e. ${{\cal P}^{*}}=\{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

[TABLE]

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

[TABLE]

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)

[TABLE]

Since

[TABLE]

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

[TABLE]

By 2) of Lemma 6.1, the function $f$ has unique maximizer $q^{*}$ (6.4). 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.3), it follows:

[TABLE]

so that with

[TABLE]

we have

[TABLE]

Since $f(q)$ is as in (6.2), it has the unique maximizer $q^{*}$ given in (6.3). On the other hand, the maximizer of $mf_{m}(q)$ is the same as the maximizer of $f_{m}(q)$, which corresponds to (6.1) 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

[TABLE]

and let $q^{*}$ be the maximizer of (6.2). 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.2. 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:

[TABLE]

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)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$,

[TABLE]

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

[TABLE]

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

1)

Suppose $\lambda=0$ and (6.6) 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

[TABLE]

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.7) holds with

[TABLE]

2)

Suppose $\lambda\in(0,1)$ and (6.7) 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.7) holds with

[TABLE]

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

[TABLE]

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.8) becomes convex. The same argument holds for (6.7). We shall present some sufficient conditions for convexity of (6.6) and (6.7) 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 [13] (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

[TABLE]

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

[TABLE]

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

[TABLE]

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.6) and (6.7) has all components strictly positive so that all maximizers of maximizer of (6.6) are interior points and the same holds for (6.7). We take $g(q)=-\ln f(q)$, where $f$ is as in (6.6) or (6.7). Thus, for any $m$, $\exp[-mg(y)]=f^{m}(q)$ so that the measure defined in the statement of theorem is

[TABLE]

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.8) is

[TABLE]

where

[TABLE]

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

[TABLE]

we obtain the Hessian

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Therefore, it is not hard to see that when

[TABLE]

then the Hessian is always positive definite, i.e. the function $g$ in (6.9) 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.10) 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.10) 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.7)

[TABLE]

and so the Hessian is

[TABLE]

where

[TABLE]

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

[TABLE]

Thus, when the following inequality holds

[TABLE]

then the Hessian is always positive definite and the function (6.11) strictly convex implying that the minimum is unique. If the Hessian is not always unique, but (6.12) holds for minimums, then Theorem 6.2 applies. Again, when $\alpha=\beta$ and $\phi(1,2)=\phi(2,1)$, then under (6.13) 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.13) holds. It means the minimum is unique, $z_{1}=z_{2}$, and one can verify that

[TABLE]

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.13) 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.13) fails. It turns out that now the function is not convex and there are two minima:

[TABLE]

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.12) holds, so that by Theorem 6.2,

[TABLE]

The function

[TABLE]

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.7) holds, with

[TABLE]

6.3 An application: exact calculation of stationary distributions for epistatic fitness

Exchangeable sampling with a product of Dirichlet priors is an abstract model of sexual reproduction, with the constraint that the population is in linkage equilibrium. For a considerable class of fitness functions it is possible to compute the stationary distribution exactly, so that the model is potentially a tool for the investigation of the properties of genetic architectures. This line of research is beyond the scope of this paper, but we give one example of computation of the stationary distribution for a non-trivial fitness function in a multi-locus system.

Consider a the case where $\mathcal{X}=\{1,2\}^{L}$, and the fitness of a genome $x=(x^{1},\ldots,x^{L})$ is defined as

[TABLE]

In other words, if a genome $(x^{1},\ldots,x^{L})$ is ‘perfect’, consisting only of the preferred alleles (represented as 1s), then it has a fitness advantage of $s$, but if even a single element is imperfect, then there is no fitness advantage. It seems plausible that there are that there are combinations of well-separated alleles in biology that have a fitness advantage only if all alleles in the set are of the correct type; this is a simple model of such a situation.

To keep the notation simple, we assume that all independent components of breeding process $\xi^{l}$, $l=1,\ldots,L$ have the same Dirichlet prior with parameter $\alpha=(\alpha_{1},\ldots,\alpha_{K})$. We shall denote the law of the breeding process by $P_{\xi}(\cdot|L,\alpha)$. When $(\xi_{1},\ldots,\xi_{n})\sim P_{\xi}(\cdot|L,\alpha)$, then the probability that the number of prefect genomes (vectors $\xi_{i}$ consisting of ones, only) equals to $n_{c}$ is

[TABLE]

In words, $P_{\xi}(n_{c},n-n_{c};L,\alpha)$ denotes the probability of $n_{c}$ perfect genomes and $n-n_{c}$ imperfect genomes in a population of $n$ genomes, each of length $L$, sampled from a product of $L$ Polya urn processes, each with prior parameter $\alpha$.

We calculate $P_{\xi}(n_{c},n-n_{c};L,\alpha)$ using recursion on $L$. When $L=1$, $n_{c}$ is simply the number of 1s in the population of $n$, so

[TABLE]

It is not hard to show, using independence of $\xi^{1},\ldots,\xi^{L}$ and exchangeability, that

[TABLE]

Finally, the stationary distribution of the number of perfect genomes in a population of $n$ is:

[TABLE]

where the normalising constant is

[TABLE]

Thus in our model the stationary distribution of the number of perfect genomes is simple to derive without making any approximations. In general, because the processes $\xi^{1},\ldots,\xi^{L}$ are independent, and as a result of the factorisation of the stationary distribution in equation 1.2, the stationary distribution for a fitness function $w$ can be derived exactly using the sum-product algorithm if $w$ can be represented as a tractable factor graph [17] over the $L$ urn distributions. This introduces a new method for studying the effect of epistatic fitness on population structure.

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.3) 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 average proportion of the first genotype has increased from 0.3 (according to the prior) to 0.9. Figure 1 illustrates the convergence.

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

The solution of (6.4) is

[TABLE]

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$.

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

[TABLE]

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

[TABLE]

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]$:

[TABLE]

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.

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

implementable computational 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, [19], but our techniques allow a proof of this effect.

Finally there is a more general conclusion from our analysis. For many years, since [15] and [12], 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

[TABLE]

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)$

[TABLE]

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.8) 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$,

[TABLE]

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

[TABLE]

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

[TABLE]

provided $n$ is big enough. Thus,

[TABLE]

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

[TABLE]

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

[TABLE]

we define Lagrangian

[TABLE]

(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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

2)

To find

[TABLE]

we define Lagrangian

[TABLE]

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

[TABLE]

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

[TABLE]

After rewriting (8.6), we obtain

[TABLE]

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

Research by Jüri Lember is supported by Estonian Institutional research funding IUT34-5 and PRG 865. Research by Chris Watkins is supported by grant number (FQXi Grant number FQXi-RFP-IPW-1913) from the Foundational Questions Institute and Fetzer Franklin Fund, a donor advised fund of Silicon Valley Community Foundation.