The Rescaled Pólya Urn and the Wright—Fisher Process with Mutation

In recent papers the authors introduce, study and apply a variant of the Eggenberger—Pólya urn, called the "rescaled" Pólya urn, which, for a suitable choice of the model parameters, exhibits a reinforcement mechanism mainly based on the last observations, a random persistent fluctuation of the predictive mean and the almost sure convergence of the empirical mean to a deterministic limit. In this work, motivated by some empirical evidence, we show that the multidimensional Wright—Fisher diffusion with mutation can be obtained as a suitable limit of the predictive means associated to a family of rescaled Pólya urns.

Keywords: Pólya urn; predictive mean; urn model; Wright—Fisher diffusion

The well-known standard Eggenberger—Pólya urn [[1]] works as follows. An urn initially contains $$N_{0,i}$$ balls of color i, for $$i=1,...,k$$ , and at each time-step, a ball is drawn from the urn and then it is returned into the urn together with $$\alpha >0$$ additional balls of the same color (here and in the following, the expression "number of balls" is not to be understood literally, but all the quantities are real numbers, not necessarily integers). Hence, denoting by $$N_{n,i}$$ the number of balls of color i inside the urn at time-step n, we have

$$N_{n,i}=N_{n-1,i}+\alpha {\xi }_{n,i}\mathrm{for}n\ge 1,$$

where $${\xi }_{n,i}=1$$ if the drawn ball at time-step n is of color i, and $${\xi }_{n,i}=0$$ otherwise. The parameter $$\alpha$$ tunes the reinforcement mechanism: the greater the $$\alpha$$ , the greater the dependence of $$N_{n,i}$$ on $${\sum }_{h=1}^n{\xi }_{h,i}$$ .

In [[3], [5]], the rescaled Pólya (RP) urn has been introduced, studied, generalized and applied. This model differs from the original one by the introduction of a parameter $$\beta$$ such that

$$\begin{array}{cccc}\hfill N_{n,i}& =b_i+B_{n,i}\hfill & \hfill & \mathrm{with}\hfill \\ \hfill B_{n+1,i}& =\beta B_{n,i}+\alpha {\xi }_{n+1,i}\hfill & \hfill & n\ge 0.\hfill \end{array}$$

Therefore, at time-step 0, the urn contains $$b_i+B_{0,i}>0$$ balls of color i and the parameters $$\alpha >0$$ and $$\beta \ge 0$$ regulate the reinforcement mechanism. More precisely, the term $$\beta B_{n,i}$$ connects $$N_{n+1,i}$$ to the "configuration" at time-step n by means of the "scaling" parameter $$\beta$$ , and the term $$\alpha {\xi }_{n+1,i}$$ connects $$N_{n+1,i}$$ to the outcome of the drawing at time-step $$n+1$$ by means of the parameter $$\alpha$$ . The case $$\beta =1$$ corresponds to the standard Eggenberger—Pólya urn with an initial number $$N_{0,i}=b_i+B_{0,i}$$ of balls of color i. When $$\beta <1$$ , the RP urn model shows the following three characteristics:

• (i) A reinforcement mechanism mainly based on the last observations;
• (ii) A random persistent fluctuation of the predictive mean $${\psi }_{n,i}=E[{\xi }_{n+1,i}=1|,{\xi }_{h,j},0\le h\le n,1\le j\le k]$$ ;
• (iii) The almost sure convergence of the empirical mean $${\sum }_{n=1}^N{\xi }_{n,i}/N$$ to the deterministic limit $$p_i=b_i/{\sum }_{i=1}^n{b}_i$$ , and a chi-squared goodness of fit result for the long-term probability distribution $$\{p_1,...,p_k\}$$ .

Regarding point (iii), we specifically have that the chi-squared statistics

$${\chi }^2=N\sum _{i=1}^k\frac{{(O_i/N-p_i)}^2}{p_i},$$

where N is the sample size and $$O_i={\sum }_{n=1}^N{\xi }_{n,i}$$ the number of sampled observations equal to i, is asymptotically distributed as $${\chi }^2(k-1)\lambda$$ , with $$\lambda >1$$ . Therefore, the presence of correlation among observations attenuates the effect of N, which multiplies the chi-squared distance between the observed frequencies and the expected probabilities. This is a key feature for statistical applications in the framework of a "big sample", where a small value of the chi-squared distance might be significant, and hence a correction related to the correlation between observations is required. In [[3], [5]], a possible application in the context of clustered data was described, with independence between clusters and correlation due to a reinforcement mechanism inside each cluster.

In [[4]], the RP urn was applied as a good model for the evolution of the sentiment associated with Twitter posts. Precisely, we analyzed three data sets: (i) the "COVID-19 epidemic" data set covers the period from 21 February to 20 April to 2020 and includes tweets in Italian about the COVID-19 epidemic; (ii) the "Migration debate" data set refers to the period from 23 January to 22 February 2019 and the collected posts are related to the Italian debate on migration; (iii) the "10 days of traffic" data set collects the entire traffic of posts in Italian in the period from 1 September to 10 September 2019. For every post, the relative sentiment, that is, the positive or negative connotation of the text, was computed using the polyglot python module developed in [[6]], which provides a numerical value $$v\in [-1,1]$$ for the sentiment of a post (for a survey on sentiment analysis, also known as opinion mining, we refer to [[7]] and references therein). We fixed a threshold T so that a tweet with $$v>T$$ was classified as a tweet with a positive sentiment and one with $$v<-T$$ was classified as a tweet with a negative sentiment. Tweets with a value $$v\in [-T,T]$$ were discarded. We took the following different values for T: $$T=0$$ , $$T=0.35$$ and $$T=0.5$$ . We applied the RP urn model, ordering the tweets according to their creation time and taking each tweet with a positive/negative classification as an extraction in the urn model. More specifically, we applied the RP model with $$k=2$$ : the time series of the tweets represents the time series of the extractions from the urn, that is, the random variables $${\xi }_{n,1}$$ . The event $$\{{\xi }_{n,1}=1\}$$ means that tweet n exhibits a positive sentiment, while $$\{{\xi }_{n,1}=0\}$$ means that tweet n exhibits a negative sentiment. For all the considered data sets, the estimated values of $$\beta$$ were strictly smaller than 1, but very near to 1 (details about the parameters estimation can be found in [[4]]). Note that the RP urn dynamics with such a value for $$\beta$$ cannot be approximated by the standard Pólya urn ( $$\beta =1$$ ), because one would lose the fluctuations of the predictive means and the possibility of touching the barriers $$\{0,1\}$$ . In this work, we show that the law of such an RP urn process can be approximated by a Wright—Fisher diffusion with mutation. More precisely, we prove that the multidimensional Wright—Fisher diffusion with mutation can be obtained as a suitable limit of the predictive means associated with a family of RP urns with $$\beta \in$$ [0, 1), $$\beta \to 1$$ . As an example, in Figure 1, for the data set "COVID-19 epidemic", we show the plot of the process $${({\psi }_{n,1})}_n$$ , reconstructed from the data (details about the reconstruction process can be found in [[4]]) and rescaled in time as $$t=n{(1-\beta )}^2$$ , the plot of a simulated (by the Euler–Maruyama method) trajectory of the Wright—Fisher process, the plot of the approximation of this trajectory by means of the RP urn and the approximation of the data process by means of the standard Pólya urn.

The Wright—Fisher (WF) class of diffusion processes models the evolution of the relative frequency of a genetic variant, or allele, in a large randomly mating population with a finite number k of genetic variants. When $$k=2$$ , the WF diffusion obeys the one-dimensional stochastic differential equation

(1) $$d{X}_t=F(X_t)dt+\sqrt{X_t(1-X_t)}d{W}_t,X_0=x_0,t\in [0,T].$$

The drift coefficient, $$F:[0,1]\to R$$ , can include a variety of evolutionary forces such as mutation and selection. For example, $$F(x)=p_1-(p_1+p_2)x=p_1(1-x)-p_{2}x$$ describes a process with recurrent mutation between the two alleles, governed by the mutation rates $$p_1>0$$ and $$p_2>0$$ . The drift vanishes when $$x=p_1/(p_1+p_2)$$ which is an attracting point for the dynamics. Equation (1) can be generalized to the case $$k>2$$ . The WF diffusion processes are widely employed in Bayesian statistics, as models for time-evolving priors [[8], [10]] and as a discrete-time finite-population construction method of the two-parameter Poisson–Dirichlet diffusion [[12]]. They have been applied in genetics [[13], [15], [17]], in biophysics [[19]], in filtering theory [[21]] and in finance [[23]].

The benefit coming from the proven limit result is twofold. First, the known properties of the WF process can give a description of the RP urn when the parameter $$\beta$$ is strictly smaller than one, but very near to one. Second, the given result might furnish the theoretical base for a new simulation method of the WF process. Indeed, the simulation from Equation (1) is highly nontrivial because there is no known closed form expression for the transition function of the diffusion, even in the simple case with null drift [[25]].

The rest of the paper is organized as follows. In Section 2, we set up our notation and we formally define the RP urn model. Section 3 provides the main result of this work, that is, the convergence result of a suitable family of predictive means associated with RP urns with $$\beta \to 1$$ . In Section 4, employing the boundary classification of the WF diffusion with mutation and connecting it to the parameters of the RP urn model, we introduce an RP urn with a value of $$\beta$$ very near to 1 the notion of recessive subsets of colors and the notion of dominant color. These two concepts are related to the possibility of reaching the barriers 0 and 1 by the predictive means of the urn process. Finally, Section 5 summarizes the work and concludes it.

For a vector $$\mathit{x}={(x_1,...,x_k)}^{\top }\in {\mathbb{R}}^k$$ , we set $$|\mathit{x}|={\sum }_{i=1}^k|x_i|$$ and $${\parallel \mathit{x}\parallel }^2={\mathit{x}}^{\top }\mathit{x}={\sum }_{i=1}^k{|x_i|}^2$$ . Moreover we denote by $$\mathbf{1}$$ and $$\mathbf{0}$$ the vectors with all the components equal to 1 and equal to 0, respectively.

Let $$\alpha >0$$ and $$\beta \ge 0$$ . At time-step 0, the urn contains $$b_i+B_{0,i}>0$$ distinct balls of color i, with $$i=1,...,k$$ . We set $$\mathit{b}={(b_1,...,b_k)}^{\top }$$ and $${\mathit{B}}_{\mathbf{0}}={(B_{0,1},...,B_{0,k})}^{\top }$$ . We suppose $$b=|\mathit{b}|>0$$ and we set $$\mathit{p}=\frac{\mathit{b}}{b}$$ . At each time-step $$(n+1)\ge 1$$ , a ball is drawn at random from the urn and we define the random vector $${\mathbf{\xi }}_{\mathit{n}+\mathbf{1}}={({\xi }_{n+1,1},...,{\xi }_{n+1,k})}^{\top }$$ as

$${\xi }_{n+1,i}=\begin{array}{cc}1\hfill & \mathrm{when}\mathrm{the}\mathrm{drawn}\mathrm{ball}\mathrm{at}\mathrm{time-step}n+1\mathrm{is}\mathrm{of}\mathrm{color}i\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}$$

The number of balls inside the urn is updated as follows:

(2) $${\mathit{N}}_{\mathit{n}+\mathbf{1}}=\mathit{b}+{\mathit{B}}_{\mathit{n}+\mathbf{1}}\mathrm{with}{\mathit{B}}_{\mathit{n}+\mathbf{1}}=\beta {\mathit{B}}_{\mathit{n}}+\alpha {\mathbf{\xi }}_{\mathit{n}+\mathbf{1}},$$

which gives

(3) $${\mathit{B}}_{\mathit{n}}={\beta }^n{\mathit{B}}_{\mathbf{0}}+\alpha {\beta }^n\sum _{h=1}^n{\beta }^{-h}{\mathbf{\xi }}_{\mathit{h}}.$$

Similarly, from the equality

$$|{\mathit{B}}_{\mathit{n}+\mathbf{1}}|=\beta |{\mathit{B}}_{\mathit{n}}|+\alpha ,$$

we get, using $${\sum }_{h=0}^{n-1}{x}^h=(1-x^n)/(1-x)$$ ,

(4) $$|{\mathit{B}}_{\mathit{n}}|={\beta }^n|{\mathit{B}}_{\mathbf{0}}|+\alpha \sum _{h=1}^n{\beta }^{n-h}={\beta }^n|{\mathit{B}}_{\mathbf{0}}|-\frac{\alpha }{1-\beta }+\frac{\alpha }{1-\beta }.$$

Setting $$r_n^{*}=|{\mathit{N}}_{\mathit{n}}|=b+|{\mathit{B}}_{\mathit{n}}|$$ , that is the total number of balls inside the urn at time-step n, we get the relations

(5) $$r_{n+1}^{*}=r_n^{*}+(\beta -1)|{\mathit{B}}_{\mathit{n}}|+\alpha$$

and

(6) $$r_n^{*}=b+\frac{\alpha }{1-\beta }+{\beta }^n|{\mathit{B}}_{\mathbf{0}}|-\frac{\alpha }{1-\beta }.$$

Denoting by $${\mathcal{F}}_0$$ the trivial $$\sigma$$ -field and setting $${\mathcal{F}}_n=\sigma ({\mathbf{\xi }}_{\mathbf{1}},...,{\mathbf{\xi }}_{\mathit{n}})$$ for $$n\ge 1$$ , the conditional probabilities $${\mathbf{\psi }}_{\mathit{n}}={({\psi }_{n,1},...,{\psi }_{n,k})}^{\top }$$ of the extraction process, also called predictive means, are

(7) $${\mathbf{\psi }}_{\mathit{n}}=E[{\mathbf{\xi }}_{\mathit{n}+\mathbf{1}}|{\mathcal{F}}_n]=\frac{{\mathit{N}}_{\mathit{n}}}{|{\mathit{N}}_{\mathit{n}}|}=\frac{\mathit{b}+{\mathit{B}}_{\mathit{n}}}{r_n^{*}}n\ge 0$$

and, from (3) and (4), we have

(8) $${\mathbf{\psi }}_{\mathit{n}}=\frac{\mathit{b}+{\beta }^n{\mathit{B}}_{\mathbf{0}}+\alpha {\sum }_{h=1}^n{\beta }^{n-h}{\mathbf{\xi }}_{\mathit{h}}}{b+\frac{\alpha }{1-\beta }+{\beta }^n|{\mathit{B}}_{\mathbf{0}}|-\frac{\alpha }{1-\beta }}.$$

The dependence of $${\mathbf{\psi }}_{\mathit{n}}$$ on $${\mathbf{\xi }}_{\mathit{h}}$$ is regulated by the factor $$f(h,n)=\alpha {\beta }^{n-h}$$ , with $$1\le h\le n,n\ge 0$$ . In the case of the standard Eggenberger—Pólya urn (i.e., the case $$\beta =1$$ ), each observation $${\mathbf{\xi }}_{\mathit{h}}$$ has the same "weight" $$f(h,n)=\alpha$$ . Instead, when $$\beta <1$$ the factor $$f(h,n)$$ increases with h, and the main contribution is given by the most recent drawings. The case $$\beta =0$$ is an extreme case, for which $${\mathbf{\psi }}_{\mathit{n}}$$ depends only on the last drawing $${\mathbf{\xi }}_{\mathit{n}}$$ .

By means of (7), together with (2) and (5), we get

(9) $${\mathbf{\psi }}_{\mathit{n}+\mathbf{1}}-{\mathbf{\psi }}_{\mathit{n}}=-\frac{(1-\beta )}{r_{n+1}^{*}}b{\mathbf{\psi }}_{\mathit{n}}-\mathit{p}+\frac{\alpha }{r_{n+1}^{*}}{\mathbf{\xi }}_{\mathit{n}+\mathbf{1}}-{\mathbf{\psi }}_{\mathit{n}}.$$

Setting $$\Delta {\mathit{M}}_{\mathit{n}+\mathbf{1}}={\mathbf{\xi }}_{\mathit{n}+\mathbf{1}}-{\mathit{\psi }}_{\mathit{n}}$$ and letting $${ϵ}_n=b(1-\beta )/r_{n+1}^{*}$$ and $${\delta }_n=\alpha /r_{n+1}^{*}$$ , from (9) we obtain

(10) $${\mathbf{\psi }}_{\mathit{n}+\mathbf{1}}-{\mathbf{\psi }}_{\mathit{n}}=-{ϵ}_n({\mathbf{\psi }}_{\mathit{n}}-\mathit{p})+{\delta }_n\Delta {\mathit{M}}_{\mathit{n}+\mathbf{1}}.$$

Consider the RP urn with parameters $$\alpha >0$$ , $$\beta \in [0,1)$$ , $$b>0$$ and $${\mathit{B}}_{\mathbf{0}}$$ such that $$|{\mathit{B}}_{\mathbf{0}}|=r(\beta )=\alpha /(1-\beta )$$ . Consequently, the total number of balls in the urn along the time-steps is constantly equal to $$r^{*}(\beta )=b+r(\beta )$$ and if we denote by $${\mathbf{\psi }}^{(\mathbf{\beta })}={({\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta })})}_n$$ the predictive means corresponding to the fixed value $$\beta$$ , we have the dynamics

(11) $${\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta })}-{\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })}=-ϵ(\beta )({\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })}-\mathit{p})+\delta (\beta )\Delta {\mathit{M}}_{\mathit{n}}^{(\mathbf{\beta })},$$

where

(12) $$ϵ(\beta )=\frac{b{(1-\beta )}^2}{\alpha +b(1-\beta )},\delta (\beta )=\frac{\alpha (1-\beta )}{\alpha +b(1-\beta )}$$

and $$\Delta {\mathit{M}}_{\mathit{n}}^{(\mathbf{\beta })}={\mathbf{\xi }}_{\mathit{n}}^{(\mathbf{\beta })}-{\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })}$$ . Note that we have $$ϵ(\beta )\sim c\delta {(\beta )}^2$$ for $$\beta \to 1$$ , with $$c=b/\alpha >0$$ . Finally, we define $${\mathit{X}}^{(\beta )}={({\mathit{X}}_{\mathit{t}}^{(\mathbf{\beta })})}_{t\ge 0}$$ , where

(13) $${\mathit{X}}_{\mathit{t}}^{(\mathbf{\beta })}={\mathbf{\psi }}_{\lfloor \mathit{t}\mathbf{/}{(\mathbf{1}-\mathbf{\beta })}^{\mathbf{2}}\rfloor }^{(\mathbf{\beta })}\iff {\mathit{X}}_{\mathit{t}}^{(\mathbf{\beta })}={\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })},t\in [(n-1){(1-\beta )}^2,n{(1-\beta )}^2).$$

The following result holds true:

Suppose that $${\mathit{X}}_{\mathbf{0}}^{(\mathbf{\beta })}$$ weakly converges towards some process $${\mathit{X}}_{\mathbf{0}}$$ when $$\beta \to 1$$ . Then, for $$\beta \to 1$$ , the family of stochastic processes $$\{{\mathit{X}}^{(\mathbf{\beta })},\beta \in [0,1)\}$$ weakly converges towards the k-alleles Wright—Fisher diffusion $$\mathit{X}={({\mathit{X}}_{\mathit{t}})}_{t\ge 0}$$ , with type-independent mutation kernel given by $$\mathit{p}$$ and with dynamics

(14) $$d{\mathit{X}}_{\mathit{t}}=-b\frac{{\mathit{X}}_{\mathit{t}}-\mathit{p}}{\alpha }dt+\Sigma ({\mathit{X}}_{\mathit{t}})d{\mathit{W}}_{\mathit{t}},$$

with $$\Sigma ({\mathit{X}}_{\mathit{t}})\Sigma {({\mathit{X}}_{\mathit{t}})}^{\top }=(\mathrm{diag}({\mathit{X}}_{\mathit{t}})-{\mathit{X}}_{\mathit{t}}{{\mathit{X}}_{\mathit{t}}}^{\top })$$ and $${\mathbf{1}}^{\top }\Sigma ({\mathit{X}}_{\mathit{t}})={\mathbf{0}}^{\top }$$ , that is,

(15) $$\begin{array}{cc}\hfill & \Sigma {({\mathit{X}}_{\mathit{t}})}_{ij}=\begin{array}{cc}0\hfill & if{X}_{t,i}{X}_{t,j}=0ori<j\hfill \\ \sqrt{X_{t,i}\frac{{\sum }_{l=i+1}^k{X}_{t,l}}{{\sum }_{l=i}^k{X}_{t,l}}}\hfill & ifi=jand{X}_{t,i}{X}_{t,j}\ne 0\hfill \\ -X_{t,i}\sqrt{\frac{X_{t,j}}{{\sum }_{l=j}^k{X}_{t,l}{\sum }_{l=j+1}^k{X}_{t,l}}}\hfill & ifi>jand{X}_{t,i}{X}_{t,j}\ne 0.\hfill \end{array}\hfill \end{array}$$

Fix a sequence $$({\beta }_n)$$ , with $${\beta }_n\in [0,1)$$ and $${\beta }_n\to 1$$ . The sequence of processes $$\{{\mathit{X}}^{({\mathbf{\beta }}_{\mathit{n}})},n\in \mathbb{N}\}$$ is bounded, hence we have to prove the tightness of the sequence in the space $$D^k[0,\infty )$$ of right-continuous functions with the usual Skorohod topology, and the characterization of the law of the unique limit process.

For any $$f\in C_b^2$$ , define

(16) $$\begin{array}{cc}\hfill {\gamma }_n^{(\beta ,f)}(\mathit{x})& ={\widehat{A}}^{(\beta )}f((n-1){(1-\beta )}^2)(\mathit{x})\hfill \\ \hfill & =E[\frac{f({\mathit{X}}_{\mathit{n}{(\mathbf{1}-\mathbf{\beta })}^{\mathbf{2}}}^{(\mathbf{\beta })})-f({\mathit{X}}_{(\mathit{n}-\mathbf{1}){(\mathbf{1}-\mathbf{\beta })}^{\mathbf{2}}}^{(\mathbf{\beta })})}{{(1-\beta )}^2}|{\mathit{X}}_{(\mathit{n}-\mathbf{1}){(\mathbf{1}-\mathbf{\beta })}^{\mathbf{2}}}^{(\mathbf{\beta })}=\mathit{x}]\hfill \\ \hfill & =E[\frac{f({\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta })})-f({\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })})}{{(1-\beta )}^2}|{\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })}=\mathit{x}]\hfill \\ \hfill & \stackrel{\mathrm{by}{\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta })}-{\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })}=-ϵ(\beta ){\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta })}-\mathit{p}+\delta (\beta )\Delta {\mathit{M}}_{\mathit{n}}^{(\mathbf{\beta })}}{=}\frac{1}{{(1-\beta )}^2}(E[f(\mathit{x})+\sum i\frac{\partial f}{\partial x_i}(\mathit{x})(-ϵ(\beta )(x_i-p_i)+\delta (\beta )\Delta {M_{n,i}}^{(\beta )})\hfill \\ \hfill & +\frac{1}{2}\delta {(\beta )}^2\sum ij\frac{{\partial }^{2}f}{\partial x_i\partial x_j}(\mathit{x})\Delta M_{n,i}^{(\beta )}\Delta M_{n,j}^{(\beta )}+O({(1-\beta )}^3)|{\mathcal{F}}_{n-1}]-f(\mathit{x}))\hfill \\ & =-\frac{b}{\alpha +b(1-\beta )}\sum i\frac{\partial f}{\partial x_i}(\mathit{x})(x_i-p_i)+\frac{1}{2}\frac{{\alpha }^2}{{(\alpha +b(1-\beta ))}^2}\sum ij\frac{{\partial }^{2}f}{\partial x_i\partial x_j}(\mathit{x})(x_i{}_{i=j}-x_i{x}_j)\hfill \\ \hfill & +O(1-\beta ).\hfill \end{array}$$

We note that, for any $$f\in C_b^2$$ , the partial derivatives in (16) are uniformly bounded, as $$\mathit{x}$$ belongs to the compact simplex $$S=\{x_i\ge 0,{\sum }_i{x}_i=1\}$$ . The family $$\{{\gamma }_n^{(\beta ,f)}(\mathit{x}),n\in \mathbb{N},\beta <1,\mathit{x}\in S\}$$ is then uniformly integrable. Thus, as a consequence of [[26]] (Theorem 4) (or [[27]] (ch. 7.4.3, Theorem 4.3, p. 236)), we have that the sequence of processes $$\{{\mathit{X}}^{({\beta }_{\mathit{n}})},n\in \mathbb{N}\}$$ is tight in the space of right-continuous functions with the usual Skorohod topology. Since, for any n and t, $${\mathit{X}}_{\mathit{t}}^{({\mathit{\beta }}_{\mathit{n}})}\in S$$ , then $${\mathbf{1}}^{\top }\Sigma ({\mathit{X}}_{\mathit{t}})={\mathbf{0}}^{\top }$$ . Moreover, the generator of the limit process is determined by the limit

$$\begin{array}{cc}\hfill Af(t)(\mathit{x})& =limn\to \infty {\gamma }_{\lfloor t/{(1-\beta )}^2\rfloor }^{({\beta }_n,f)}(\mathit{x})\hfill \\ & =-\frac{b}{\alpha }\sum i\frac{\partial f}{\partial x_i}(\mathit{x})(x_i-p_i)+\frac{1}{2}\sum ij\frac{{\partial }^{2}f}{\partial x_i\partial x_j}(\mathit{x})(x_i{}_{i=j}-x_i{x}_j).\hfill \end{array}$$

Hence, the weak limit of the sequence of the bounded processes $${\mathit{X}}^{({\mathit{\beta }}_{\mathit{n}})}$$ is the diffusion process

$$d{\mathit{X}}_{\mathit{t}}=-b\frac{{\mathit{X}}_{\mathit{t}}-\mathit{p}}{\alpha }dt+\Sigma ({\mathit{X}}_{\mathit{t}})d{\mathit{W}}_{\mathit{t}},\Sigma ({\mathit{X}}_{\mathit{t}})\Sigma {({\mathit{X}}_{\mathit{t}})}^{\top }=(\mathrm{diag}({\mathit{X}}_{\mathit{t}})-{\mathit{X}}_{\mathit{t}}{{\mathit{X}}_{\mathit{t}}}^{\top }).$$

The expression (15) follows from [[28]] (Corollary 3). □

Since the simplex has dimension $$k-1$$ with respect to the Lebesgue measure, it is convenient to change the notations. Let $$T^{k-1}$$ be the $$k-1$$ -dimensional simplex defined by

$$T^{k-1}:=\{\mathit{y}\in {\mathbb{R}}^{k-1}:y_1\ge 0,...,y_{k-1}\ge 0,1-y_1-y_2-...-y_{k-1}\ge 0\},$$

where, with the old definition, we have $$x_i=y_i,i<k$$ and $$x_k:=1-y_1-y_2-...-y_{k-1}$$ . Obviously, there is a one-to-one natural correspondence between $$T^{k-1}$$ and the simplex $$\{\mathit{x}\in {\mathbb{R}}^k:x_1\ge 0,...,x_k\ge 0,{\sum }_i{x}_i=1\}$$ defined by

$$\mathit{y}=(y_1,...,y_{k-1})⟷(y_1,...,y_{k-1},1-y_1-y_2-...-y_{k-1})=(x_1,...,x_{k-1},x_k)=\mathit{x}.$$

The Markov diffusion process $${\mathit{X}}_{\mathit{t}}$$ in (14) may be redefined as $${\mathit{Y}}_{\mathit{t}}=(X_{t,1},...,X_{t,k-1})$$ on $$\mathit{y}\in T^{k-1}$$ with the corresponding generator

(17) $$Lf(\mathit{y})=-\frac{b}{\alpha }\sum _{i=1}^{k-1}\frac{\partial f}{\partial y_i}(\mathit{y})(y_i-p_i)+\frac{1}{2}\sum _{i,j=1}^{k-1}\frac{{\partial }^{2}f}{\partial y_i\partial y_j}(\mathit{y})(y_i{}_{i=j}-y_i{y}_j).$$

The Kolmogorov forward equation for the density $$p(\mathit{y},t)$$ of the limiting process $${\mathit{Y}}_t$$ is

(18) $$\begin{array}{c}\frac{\partial }{\partial t}p(\mathit{y},t)=\frac{1}{2}(\frac{b}{\alpha }\sum _{i=1}^{k-1}\frac{\partial }{\partial y_i}(p(\mathit{y},t)(y_i-p_i))\hfill \\ +\sum _{i=1}^{k-1}\frac{{\partial }^2}{\partial y_i^2}(y_i(1-y_i)p(\mathit{y},t))-2\sum 1\le i<j\le k-1\frac{{\partial }^2}{\partial y_i\partial y_j}(y_i{y}_{j}p(\mathit{y},t))).\hfill \end{array}$$

Therefore, it is not hard to show that the limit invariant ergodic distribution is

(19) $$p(\mathit{y})=\frac{1}{B(2\frac{b}{\alpha }\mathit{p})}{(1-y_1-...-y_{k-1})}^{\frac{2b(1-p_1-...-p_{k-1})}{\alpha }-1}\prod _{i=1}^{k-1}{y}_i^{\frac{2b{p}_i}{\alpha }-1},$$

because it satisfies (18) (see also [[29]]). The above distribution is the Dirichlet distribution $$\mathit{Dir}2\frac{b}{\alpha }\mathit{p}$$ as a function of $$\mathit{x}=(\mathit{y},1-y_1-...-y_{k-1})$$ .

The transition density $$p({\mathit{y}}_{\mathbf{0}},\mathit{y};t)$$ is defined by

$$P({\mathit{Y}}_{\mathit{t}}\in S|{\mathit{Y}}_{\mathbf{0}}={\mathit{y}}_{\mathbf{0}})={\int }_{S\cap T^{k-1}}p({\mathit{y}}_{\mathbf{0}},\mathit{y};t)d\mathit{y}$$

and it can be represented in terms of series of orthogonal polynomials [[30]] as shown in [[31]]. Moreover, we refer to [[9],[32],[33]] for the explicit form of the reproducing kernel orthogonal polynomials.

Let $$\mathit{J}=\{J_1,...,J_{k_J}\}$$ be a partition of $$\{1,...,k\}$$ , in that $$J_l\ne Ø$$ , $$J_{i_1}\cap J_{i_2}=Ø$$ , and $${\cup }_{l=1}^{k_J}=\{1,...,k\}$$ . Here $$k_j$$ denotes the cardinality of $$\mathit{J}$$ . Define the $$k_J$$ -dimensional objects $${({\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta },\mathit{J})})}_n$$ , $${({\mathbf{\xi }}_{\mathit{n}}^{(\mathbf{\beta },\mathit{J})})}_n$$ and $${\mathit{p}}^{(\mathit{J})}$$ as

$$\begin{array}{cc}\hfill {\psi }_{n,i}^{(\beta ,\mathit{J})}& =\sum l\in J_i{\psi }_{n,l}^{(\beta )}\hfill \\ \hfill {\xi }_{n,i}^{(\beta ,\mathit{J})}& =\sum l\in J_l{\xi }_{n,l}^{(\epsilon )}\hfill \\ \hfill p_i^{(\mathit{J})}& =\sum l\in J_i{p}_l\hfill \end{array}\mathrm{for}i=1,...,k_J,$$

and $${\mathit{X}}_{\mathit{t}}^{(\mathbf{\beta },\mathit{J})}={\mathbf{\psi }}_{\lfloor \mathit{t}\mathbf{/}{(\mathbf{1}-\mathbf{\beta })}^{\mathbf{2}}\rfloor }^{(\mathbf{\beta },\mathit{J})}$$ . With these definitions, from (11), we immediately get that $${({\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta },\mathit{J})})}_n$$ is a $$k_J$$ -dimensional RP urn following the dynamics

(20) $${\mathbf{\psi }}_{\mathit{n}}^{(\mathbf{\beta },\mathit{J})}-{\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta },\mathit{J})}=-ϵ(\beta )({\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta },\mathit{J})}-{\mathit{p}}^{(\mathit{J})})+\delta (\beta )({\mathbf{\xi }}_{\mathit{n}}^{(\mathbf{\beta },\mathit{J})}-{\mathbf{\psi }}_{\mathit{n}-\mathbf{1}}^{(\mathbf{\beta },\mathit{J})})$$

and that Theorem 1 holds for $${\mathit{X}}_{\mathit{t}}^{(\mathbf{\beta },\mathit{J})}$$ . Consequently, the convergence to the Wright—Fisher diffusion still holds if we group together some components of the process. For instance, when we consider two groups of components, we have the following result:

Let $$\mathit{J}=\{J,J^c\}$$ with $$J\ne Ø$$ , $$J^c\ne Ø$$ . Under the hypothesis of Theorem 1, each component of the sequence of processes $${\mathit{X}}_{\mathit{t}}^{(\mathbf{\beta },\mathit{J})}$$ converges, for $$\beta \to 1$$ , to the one-dimensional diffusion process with values in $$[0,1]$$ that satisfies the SDE

$$d{X}_{t,i}^{(\mathit{J})}=-b\frac{X_{t,i}^{(\mathit{J})}-p_i}{\alpha }dt+{(-1)}^{i+1}\sqrt{X_{t,i}^{(\mathit{J})}(1-X_{t,i}^{(\mathit{J})})}d{W}_t.$$

In addition, $$X_{t,1}^{(\mathit{J})}={\sum }_{l\in J}{X}_{t,l}$$ and $$X_{t,2}^{(\mathit{J})}={\sum }_{l\in J^c}{X}_{t,l}$$ .

Now, if we further specialize the grouping choice to $$\mathit{J}=(\{i\},\{1,...,i-1,i+1,...,k\})$$ , we get:

Under the conditions of Theorem 1 the i-th component of the sequence of processes $${\mathit{X}}^{(\mathbf{\beta })}$$ converges, for $$\beta \to 1$$ , to the one-dimensional diffusion $${(X_{t,i})}_{t\ge 0}$$ with values in $$[0,1]$$ satisfying the SDE

$$d{X}_{t,i}=-b\frac{X_{t,i}-p_i}{\alpha }dt+\sqrt{X_{t,i}(1-X_{t,i})}d{W}_t.$$

For instance, the above two results are useful in order to translate the well-known classification of the boundaries of the WF process with mutation [[34]] (p. 239, Example 8) (see also [[35]]) to the RP urn model when the parameter $$\beta$$ is strictly smaller than 1, but very near to 1. Indeed, Corollary 1 implies that $$Z_t={\sum }_{l\in J}{X}_{t,l}$$ satisfies the SDE

$$\begin{array}{cc}\hfill d{Z}_t& =-b\frac{Z_t-{\sum }_{l\in J}{p}_l}{\alpha }dt+\sqrt{Z_t(1-Z_t)}d{W}_t\hfill \\ & =-\frac{b}{\alpha }(1-\sum l\in J{p}_l)Z_t+\frac{b}{\alpha }\sum l\in J{p}_l(1-Z_t)dt+\sqrt{Z_t(1-Z_t)}d{W}_t.\hfill \end{array}$$

Setting $$a_0=\frac{b}{\alpha }{\sum }_{l\in J}{p}_l$$ and $$a_1=\frac{b}{\alpha }-a_0$$ and noting that $${\cap }_{i\in J}\{X_{t,i}=0\}=\{Z_t=0\}$$ , we obtain:

• (1) $$a_0<1/2$$ , i.e., $${\sum }_{l\in J}{p}_l<\frac{\alpha }{2b}$$ , if and only if $$P(\exists t:{\cap }_{i\in J}\{X_{t,i}=0\})=1$$ ;
• (2) $$a_0\ge 1/2$$ , i.e., $${\sum }_{l\in J}{p}_l\ge \frac{\alpha }{2b}$$ , if and only if $$P(\exists t:{\cap }_{i\in J}\{X_{t,i}=0\})=0$$ .

With the same spirit, Corollary 2 states that $$Z_t=1-X_{t,i}$$ satisfies the SDE

$$\begin{array}{cc}\hfill d{Z}_t& =-b\frac{Z_t-{\sum }_{l\ne i}{p}_l}{\alpha }dt+\sqrt{(1-Z_t)Z_t}d{W}_t\hfill \\ & =-\frac{b}{\alpha }{p}_i{Z}_t+\frac{b}{\alpha }(1-p_i)(1-Z_t)dt+\sqrt{Z_t(1-Z_t)}d{W}_t.\hfill \end{array}$$

Setting $$a_0=\frac{b}{\alpha }(1-p_i)$$ and $$a_1=\frac{b}{\alpha }-a_0$$ , we get:

• (3) $$a_0<1/2$$ , i.e., $$p_i>1-\frac{\alpha }{2b}$$ , if and only if $$P(\exists t:\{X_{t,i}=1\})=1$$ ;
• (4) $$a_0\ge 1/2$$ , i.e., $$p_i\le 1-\frac{\alpha }{2b}$$ , if and only if $$P(\exists t:\{X_{t,i}=1\})=0$$ .

Therefore, for an RP urn with $$\beta <1$$ , but very near to 1, we can give the following definition:

We call recessive a non-empty subset $$J⊊\{1,...,k\}$$ of colors such that $${\sum }_{l\in J}{p}_l<\frac{\alpha }{2b}$$ . We call dominant a color $$i\in \{1,...,k\}$$ such that $$\{1,...,k\}\backslash \left\{i\right\}$$ is recessive.

Obviously, every subset of a recessive set is recessive. Moreover, when $$\frac{\alpha }{b}>2(1-{min}_i{p}_i)$$ , every set $$J⊊\{1,...,k\}$$ is recessive. The terms "recessive" and "dominant" are justified by the fact that, recalling properties (1)–(4) of the WF process, if a set of colors is recessive, then we can observe that at some times the corresponding predictive means of the urn process are very near to zero. On the contrary, when a color is dominant, we can observe that at some times the corresponding predictive mean of the urn process is very near to one. In Figure 2, we plot the process $$({\psi }_{n,1})$$ related to the simulation of an RP urn with $$k=2$$ , $$\alpha /b=1$$ and $$p=0.75$$ , where it is possible to observe the excursions near the barrier 1.

We have proven that the multidimensional WF diffusion with mutation can be obtained as the limit of the predictive means associated with a family of RP urns with $$\beta <1$$ , $$\beta \to 1$$ . As a consequence, the known properties of the WF process can give a description of the RP urn when the parameter $$\beta$$ is strictly smaller than 1, but very near to 1. For instance, starting from the known classification of the boundaries for the WF process and connecting it to the model parameters of the RP urn, we have obtained for an RP urn with a value of $$\beta$$ very near to one, the notion of recessive subsets of colors and the notion of a dominant color. These two concepts are related to the possibility of reaching the barriers 0 and 1 by the predictive means of the urn process. Other classical problems, together with the corresponding known results for the WF process, can be found in [[31]]. These results can be used in order to give an approximated answer to the considered problems in the case of an RP urn with a value of $$\beta$$ near 1.

Graph: Figure 1 "COVID-19 epidemic" Twitter data set: the black line is the process (ψn,1)n , reconstructed from the data and rescaled in time as t=n(1−β)2 ; the red line is a simulated trajectory of the Wright—Fisher process; the orange line is the approximation of this trajectory by means of the RP urn and the blue line is the approximation of the data process by means of the standard Pólya urn. The numbers 0, 0.35 and 0.5 refer to the values chosen for the threshold T. The corresponding estimated values for 1−β are: 0.000776 (8×10−4), 0.00115 (11×10−4) and 0.00130 (13×10−4).

Graph: Figure 2 Simulation: plot of the process (ψn,1) related to the simulation of an RP urn with k=2 , α/b=1 and p=0.75.

Both authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Irene Crimaldi is partially supported by the Italian "Programma di Attività Integrata" (PAI), project "TOol for Fighting FakEs" (TOFFE) funded by the IMT School for Advanced Studies Lucca.

Not applicable.

Not applicable.

The authors declare no conflict of interest.

Both authors sincerely thank the organizers of the present special issue for their invitation to contribute and Fabio Saracco for having collected and shared with them the analyzed Twitter data sets. Giacomo Aletti is a member of the Italian group "Gruppo Nazionale per il Calcolo Scientifico" of the Italian Institute "Istituto Nazionale di Alta Matematica". Irene Crimaldi is a member of the Italian group "Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni" of the Italian Institute "Istituto Nazionale di Alta Matematica".