Adaptive Learning in Large Populations

TL;DR

This paper models adaptive learning in large populations using a PDE derived from Harley's rule, revealing that faster memory decay benefits evolution in certain 2x2 games.

Contribution

It introduces a PDE framework for stochastic adaptive learning in large populations, connecting behavioral rules to population dynamics.

Findings

Faster memory decay provides an evolutionary advantage in some games.

The PDE model captures the conservation of stimuli in behavior selection.

Analysis of 2x2 games shows conditions favoring different learning speeds.

Abstract

We consider the adaptive learning rule of Harley (1981) for behavior selection in symmetric conflict games in large populations. The rule uses organisms' past, accumulated rewards as the predictor for the future behavior, and can be traced in many life forms from bacteria to humans. We derive a partial differential equation (PDE) that describes the stochastic learning in a population of agents. The equation has simple structure of the `conservation of mass'-type equation in the space of stimuli to engage in a particular type of behavior. We analyze the solutions of the PDE model for typical 2x2 games. It is found that in games with small residual stimuli, adaptive learning rules with faster memory decay have an evolutionary advantage.

Tables

Table 1. Table 1: The incremental payoff matrix for the game. Parameters a , b , c , d , h 𝑎 𝑏 𝑐 𝑑 ℎ a,b,c,d,h are positive.

	A	B
A	(ah, ah)	(dh, ch)
B	(ch, dh)	(bh,  bh)

Equations

$$\frac{d{m}_i}{dt}=k{m}_i(W_i-\bar{W}),i=\mathrm{1..}N$$

$$\frac{S_i(t)}{{\sum }_j{S}_j(t)}.$$

$$S_{i_0}(t+1)=(1-\bar{\mu })S_{i_0}(t)+\bar{\mu }{r}_{i_0}+P_{i_0}(t),$$

$$S_k(t+1)=(1-\bar{\mu })S_k(t)+\bar{\mu }{r}_k.$$

$$X^{t+\delta}_{i}{}={}\left\{\begin{array}[]{ll}(r_{1}\mu h{}+{}(1-\mu h)S^{t}_{1,i}+ah,\,r_{2}\mu h+(1-\mu h))S^{t}_{2,i}),&\,\rm(A,A)\\ (r_{1}\mu h{}+{}(1-\mu h)S^{t}_{1,i}+ch,\,r_{2}\mu h+(1-\mu h))S^{t}_{2,i}),&\,\rm(A,B)\\ (r_{1}\mu h+(1-\mu h))S^{t}_{1,i},\,r_{2}\mu h{}+{}(1-\mu h)S^{t}_{2,i}+dh),&\,\rm(B,A)\\ (r_{1}\mu h+(1-\mu h))S^{t}_{1,i},\,r_{2}\mu h{}+{}(1-\mu h)S^{t}_{2,i}+bh),&\,\rm(B,B)\end{array}\right.$$

$$\frac{\partial f}{\partial t}+\text{mbox}\mathrm{div}(u(x,t)f)=0,$$

$$u(x,t){}={}\frac{1}{s_{1}+s_{2}}\left[\begin{array}[]{l}\left(a\bar{m}(t){}+{}d(1-\bar{m}(t)\right)s_{1}+\mu(r_{1}-s_{1})(s_{1}+s_{2})\\ \left(c\bar{m}(t){}+{}b(1-\bar{m}(t)\right)s_{2}+\mu(r_{2}-s_{2})(s_{1}+s_{2})\end{array}\right],$$

$$\bar{m}(t)=\int \frac{s_1}{s_1+s_2}f(x,t)dx.$$

$$\frac{dx}{dt}=u(x,t),$$

$$\frac{d\bar{m}}{dt}=\left((a-c)\bar{m}+(d-b)(1-\bar{m})\right)\int \frac{s_1{s}_2}{(s_1+s_2{)}^3}f(x,t)dx.$$

$$\frac{d(\bar{m}-m_{\ast })}{dt}=(c-a+d-b)(\bar{m}-m_{\ast })\int \frac{s_1{s}_2}{(s_1+s_2{)}^3}f(x,t)dx.$$

$$\int \frac{s_1}{s_1+s_2}{f}_0(x)dx=m_{\ast },$$

$$\int \frac{s_1}{s_1+s_2}f(x,t)dx=m_{\ast }.$$

$$X_i^t=(s_1(t),s_2(t)),\forall i=\mathrm{1..}N,$$

$$\frac{dm}{dt}=\frac{1}{s_1(t)+s_2(t)}m(1-m)\left((a-c)m+(d-b)(1-m)\right).$$

$$X^{t+\delta}_{i}{}={}\left\{\begin{array}[]{ll}((1-\mu h)s_{1,i}+r_{1}\mu h+ah,\,(1-\mu h)s_{2,i}+r_{2}\mu h)&\rm Prob=\lambda_{i}\lambda_{j}\\ ((1-\mu h)s_{1,i}+r_{1}\mu h+dh,\,(1-\mu h)s_{2,i}+r_{2}\mu h)&\rm Prob=\lambda_{i}(1-\lambda_{j})\\ ((1-\mu h)s_{1,i}+r_{1}\mu h,\,(1-\mu h)s_{2,i}+r_{2}\mu h+ch)&\rm Prob=(1-\lambda_{i})\lambda_{j}\\ ((1-\mu h)s_{1,i}+r_{1}\mu h,\,(1-\mu h)s_{2,i}+r_{2}\mu h+bh)&\rm Prob=(1-\lambda_{i})(1-\lambda_{j})\end{array}\right.$$

$$\int\phi(\bar{x})w_{h}(\bar{x},t+\delta)\,d\bar{x}{}={}\mathbb{E}[\phi(X^{t+h})]\\ {}={}\sum_{i\not=j}(N(N-1))^{-1}\int\left(\lambda_{i}\lambda_{j}\phi(\bar{x})\Big{|}_{x_{i}{}={}((1-\mu h)s_{1,i}+r_{1}\mu h+ah,\,(1-\mu h)s_{2,i}+r_{2}\mu h)\atop\,\,\,x_{j}=((1-\mu h)s_{1,j}+r_{1}\mu h+ah,\,(1-\mu h)s_{2,j}+r_{2}\mu h)}\right.\\ \left.+\lambda_{i}(1-\lambda_{j})\phi(\bar{x})\Big{|}_{x_{i}{}={}((1-\mu h)s_{1,i}+r_{1}\mu h+dh,\,(1-\mu h)s_{2,i}+r_{2}\mu h)\atop\,\,\,x_{j}{}={}((1-\mu h)s_{1,j}+r_{1}\mu h,\,(1-\mu h)s_{2,j}+r_{2}\mu h+ch)}\right.\\ \left.+(1-\lambda_{i})\lambda_{j}\phi(\bar{x})\Big{|}_{x_{i}=((1-\mu h)s_{1,i}+r_{1}\mu h,\,(1-\mu h)s_{2,i}+r_{2}\mu h+ch)\atop\,\,\,x_{j}=((1-\mu h)s_{1,j}+r_{1}\mu h+dh,\,(1-\mu h)s_{2,j}+r_{2}\mu h)}\right.\\ \left.+(1-\lambda_{i})(1-\lambda_{j})\phi(\bar{x})\Big{|}_{x_{i}=((1-\mu h)s_{1,i}+r_{1}\mu h,\,(1-\mu h)s_{2,i}+r_{2}\mu h+bh)\atop\,\,\,x_{j}=((1-\mu h)s_{1,j}+r_{1}\mu h,\,(1-\mu h)s_{2,j}+r_{2}\mu h+bh)}\right)w_{h}(\bar{x},t)\,d\bar{x}.$$

$$\int\phi(\bar{x})w_{h}(\bar{x},t+\delta)\,d\bar{x}{}={}\int\phi(\bar{x})w_{h}(x,t)\,d\bar{x}\\ {}+{}\sum_{i\not=j}(N(N-1))^{-1}\int\left(\lambda_{i}\lambda_{j}[\phi(\bar{x})\Big{|}_{x_{i}{}={}((1-\mu h)s_{1,i}+r_{1}\mu h+ah,\,(1-\mu h)s_{2,i}+r_{2}\mu h)\atop\,\,\,x_{j}=((1-\mu h)s_{1,j}+r_{1}\mu h+ah,\,(1-\mu h)s_{2,j}+r_{2}\mu h)}-\phi(\bar{x})]\right.\\ \left.+\lambda_{i}(1-\lambda_{j})[\phi(\bar{x})\Big{|}_{x_{i}{}={}((1-\mu h)s_{1,i}+r_{1}\mu h+dh,\,(1-\mu h)s_{2,i}+r_{2}\mu h)\atop\,\,\,x_{j}{}={}((1-\mu h)s_{1,j}+r_{1}\mu h,\,(1-\mu h)s_{2,j}+r_{2}\mu h+ch)}-\phi(\bar{x})]\right.\\ \left.+(1-\lambda_{i})\lambda_{j}[\phi(\bar{x})\Big{|}_{x_{i}=((1-\mu h)s_{1,i}+r_{1}\mu h,\,(1-\mu h)s_{2,i}+r_{2}\mu h+ch)\atop\,\,\,x_{j}=((1-\mu h)s_{1,j}+r_{1}\mu h+dh,\,(1-\mu h)s_{2,j}+r_{2}\mu h)}-\phi(\bar{x})]\right.\\ \left.+(1-\lambda_{i})(1-\lambda_{j})[\phi(\bar{x})\Big{|}_{x_{i}=((1-\mu h)s_{1,i}+r_{1}\mu h,\,(1-\mu h)s_{2,i}+r_{2}\mu h+bh)\atop\,\,\,x_{j}=((1-\mu h)s_{1,j}+r_{1}\mu h,\,(1-\mu h)s_{2,j}+r_{2}\mu h+bh)}-\phi(\bar{x})]\right)w_{h}(\bar{x},t)\,d\bar{x}.$$

$$f_{h}(x,t){}={}\sum_{k}N^{-1}\int w_{h}(\bar{x})\big{|}_{x_{k}=x}\,d\bar{x}_{k},\quad x\in\mathbb{R}^{2},$$

$$g_{h}(x,y,t){}={}\sum_{i\not=j}(N(N-1))^{-1}\int w_{h}(\bar{x})\big{|}_{x_{i}=x,\,\\ x_{j}=y}\,d\bar{x}_{ij},$$

$$f_h(x,t)=\int g_h(x,y,t)dx=\int g_h(x,y,t)dy.$$

$$\int \psi (x)f_h(x,t)dx=\sum _k{N}^{-1}\int \psi (x_k)w_h(\bar{x})d\bar{x},$$

$$\int \omega (x,y)g_h(x,y,t)dxdy=\sum _{i\ne j}(N(N-1){)}^{-1}\int \omega (x_i,x_j)w_h(\bar{x})d\bar{x}.$$

$$\begin{gathered} \int \psi (x)f_h(x,t+\delta )dx=\int \psi (x)f_h(x,t)dx \\ +N^{-1}\sum _{i\ne j}(N(N-1){)}^{-1}\int ({\lambda }_i{\lambda }_j[\psi ((1-\mu h)s_{1,i}+r_1\mu h+ah,(1-\mu h)s_{2,i}+r_2\mu h)-\psi (x_i) \\ +\psi ((1-\mu h)s_{1,j}+r_1\mu h+ah,(1-\mu h)s_{2,j}+r_2\mu h)-\psi (x_j)] \\ +{\lambda }_i(1-{\lambda }_j)[\psi ((1-\mu h)s_{1,i}+r_1\mu h+dh,(1-\mu h)s_{2,i}+r_2\mu h)-\psi (x_i) \\ +\psi ((1-\mu h)s_{1,j}+r_1\mu h,(1-\mu h)s_{2,j}+r_2\mu h+ch)-\psi (x_j)] \\ +(1-{\lambda }_i){\lambda }_j[\psi ((1-\mu h)s_{1,i}+r_1\mu h,(1-\mu h)s_{2,i}+r_2\mu h+ch)-\psi (x_i) \\ +\psi ((1-\mu h)s_{1,j}+r_1\mu h+dh,(1-\mu h)s_{2,j}+r_2\mu h)-\psi (x_j)] \\ +(1-{\lambda }_i)(1-{\lambda }_j)[\psi ((1-\mu h)s_{1,i}+r_1\mu h,(1-\mu h)s_{2,i}+r_2\mu h+bh)-\psi (x_i) \\ +\psi ((1-\mu h)s_{1,j}+r_1\mu h,(1-\mu h)s_{2,j}+r_2\mu h+bh)-\psi (x_j)])w_h(\bar{x},t)d\bar{x}. \end{gathered}$$

$$\begin{gathered} \int \psi (x)f_h(x,t+\delta )dx=\int \psi (x)f_h(x,t)dx \\ +2N^{-1}\int (\lambda (x)\lambda (y)[\psi ((1-\mu h)s_1^x+r_1\mu h+ah,(1-\mu h)s_2^x+r_2\mu h)-\psi (x) \\ +\psi ((1-\mu h)s_1^y+r_1\mu h+ah,(1-\mu h)s_2^y+r_2\mu h)-\psi (y)] \\ +\lambda (x)(1-\lambda (y))[\psi ((1-\mu h)s_1^x+r_1\mu h+dh,(1-\mu h)s_2^x+r_2\mu h)-\psi (x) \\ +\psi ((1-\mu h)s_1^y+r_1\mu h,(1-\mu h)s_2^y+r_2\mu h+ch)-\psi (y)] \\ +(1-\lambda (x))\lambda (y)[\psi ((1-\mu h)s_1^x+r_1\mu h,(1-\mu h)s_2^x+r_2\mu h+ch)-\psi (x) \\ +\psi ((1-\mu h)s_1^y+r_1\mu h+dh,(1-\mu h)s_2^y+r_2\mu h)-\psi (y)] \\ +(1-\lambda (x))(1-\lambda (y))[\psi ((1-\mu h)s_1^x+r_1\mu h,(1-\mu h)s_2^x+r_2\mu h+bh)-\psi (x) \\ +\psi ((1-\mu h)s_1^y+r_1\mu h,(1-\mu h)s_2^y+r_2\mu h+bh)-\psi (y)])g_h(x,y,t)dxdy. \end{gathered}$$

$$g_h(x,y,t)=f_h(x,t)f_h(y,t).$$

$$\begin{gathered} \frac{N}{2h}\int \psi (x)(f_h(x,t+\delta )-f_h(x,t))dx= \\ 2{\bar{m}}_h(t)\int \left(\lambda (x)(a+\mu (r_1-s_1)){\partial }_{s_1}\psi (x)+\lambda (x)\mu (r_2-s_2){\partial }_{s_2}\psi (x)\right)f_h(x,t)dx \\ +2(1-{\bar{m}}_h(t))\int \left(\lambda (x)(d+\mu (r_1-s_1)){\partial }_{s_1}\psi (x)+\lambda (x)\mu (r_2-s_2){\partial }_{s_2}\psi (x)\right)f_h(x,t)dx \\ +2{\bar{m}}_h(t)\int ((1-\lambda (x))(\mu (r_1-s_1)){\partial }_{s_1}\psi (x) \\ +(1-\lambda (x))(c+\mu (r_2-s_2)){\partial }_{s_2}\psi (x))f_h(x,t)dx \\ +2(1-{\bar{m}}_h(t))\int ((1-\lambda (x))(\mu (r_1-s_1)){\partial }_{s_1}\psi (x) \\ +(1-\lambda (x))(b+\mu (r_2-s_2)){\partial }_{s_2}\psi (x))f_h(x,t)dx \end{gathered}$$

$$\begin{gathered} \frac{N}{4h}\int \psi (x)(f_h(x,t+\delta )-f_h(x,t))dx= \\ \int \left(v_1{\partial }_{s_1}\psi +v_2{\partial }_{s_2}\psi \right)f_h(x,t)dx, \end{gathered}$$

$$v_1=\lambda (x)(a{\bar{m}}_h(t)+(1-{\bar{m}}_h(t))d)+\mu (r_1-s_1),$$

$$v_2=(1-\lambda (x))(c{\bar{m}}_h(t)+(1-{\bar{m}}_h(t))b)+\mu (r_2-s_2).$$

Full text

††journal: Journal Mathematical Biology

1 Introduction

The seminal paper Maynard Smith & Price (1973) introduced the concepts of game theory into the study of animal behavior to explain the evolution of behavioral traits. To this end, a notion of evolutionarily stable strategy (ESS) was developed to describe stable outcomes of natural selection, by defining it as being uninvadable by pure or mixed strategies players. A dynamical process leading to an ESS can be formalized by the following model. Consider a situation where each individual in a population consistently uses one of the available behavioral traits $t_{1},..,t_{k}$ to interact with other members and no mutations take place. Assume the choice of the behavior (action) is inherited, and the population consists of groups that use a particular behavior. Over a long period of time involving large number of interactions, the average fitness per game for individuals using $t_{i},$ denoted by $W_{i},$ is compared to the population average fitness $\bar{W}$, and the individuals in $t_{i}$-group are reproduced at rate proportional to $W_{i}/\bar{W}.$ If the changes in the frequencies of $t_{i}$ are approximately continuous Taylor & Joniker (1978), Zeeman (1981) derived the replicator dynamics equations for the frequencies $m_{i}$ of individuals acting according to $t_{i}:$

[TABLE]

and showed that ESS’s are the asymptotically stable fixed points of this system of equations.

Considered from the point of view of game theory, an evolutionarily stable strategy is a refinement of a Nash equilibrium, which describes an optimal choice of actions in games. This way, natural selection is a mechanism for implementing a rational decision-making in the evolution of species. There is another way by which organisms, even without complex cognition, can discover optimal actions. It can be achieved through their ability to regulate behaviors depending on the experience, in particular, through the tendency to repeat positive and avoid negative experiences. This is known as the law of effect, first formulated by Thorndike (1989), and generally accepted as one of the main paradigms of animal behavior, see for example, Ferster & Skinner (1957), Herrnstein (1961), Herrnstein (1970), Catania (1963), Chung & Herrnstein (1967), Domjan & Burkhard (1986).

The law of effect is the basis for the reinforcement learning models. They were introduced by Bush & Mosteller (1955), and since then have been applied to problems in such diverse areas as biology, economics, and engineering. Some representative examples of the extensive literature on this subject can be found in Harley (1981), Cross (1983), Roth & Erev (1995), Erev & Roth (1998), Sutton & Barto (1998), Sandholdm (2010), Nax & Perc (2015).

Interestingly enough, Börgers & Sarin (1997) demonstrated that models of learning also lead to the replicator dynamics equations similar to (1). For that, they considered repeated plays of a 2x2 game between two agents who adjust their probabilities for actions according to a reinforcement learning model of Cross (1983), and derived the replicator dynamics equation in the limit of small payoffs. More general situations, both in terms of the type of the games, the kind of reinforcement models and their relation to replicator-like equations are discussed in a monograph by Fudenberg & Levine (1998), Rustichini (1999), Fudenberg & Takahashi (2011), and Mertikopoulos & Sandholm (2016).

The law of effect can be expressed in many different ways, depending on which decision-making facilities are reinforced and on the specific rules of reinforcement. In the model of Cross (1983), it is the probability to play particular action that undergoes reinforcement. In an alternative model motivated by bio-chemical processes in neural circuits, Harley (1981) proposed the relative payoff sum (RPS) algorithm as a decision making mechanism. The RPS learning rule assumes the ability of organism to maintain a record of cumulative rewards from previous experiences, which at epoch $t$ is given by vector $S(t)=(S_{1}(t),..,S_{k}(t)).$ $S(t)$ is the predictor for future behavior and can be interpreted as a vector of “motivations” or “stimuli” to engage in a corresponding action. From the current stimuli an agent computes the probability to play $t_{i}:$

[TABLE]

If action $t_{i_{0}}$ is chosen and it brings payoff $P_{i_{0}}(t),$ which is assumed to be non-negative, then stimuli are updated according to the rule

[TABLE]

and for $k\not=i_{0},$

[TABLE]

Positive parameter $\bar{\mu}$ expresses memory effect: payoff $P_{i}(t-n)$ from $n$ previous plays will appear with the weight $(1-\bar{\mu})^{n+1}$ in the expression for $S_{i}(t+1).$ Parameter $r_{i}$ is some default level (residual) of stimuli $S_{i}.$ For example, the residuals might represent the genetic preference for the action. A similar learning model was introduced by Roth & Erev (1995).

The goal of this paper is to address the question of behavior of RPS learning agents in large populations, where agents are randomly matched in pairwise encounters, i.e. learning in heterogeneous populations. Unlike the evolutionary dynamics framework, RPS learning agents can not be identified with some particular strategy they use all the time. The RPS rule will in general prescribe new strategy every time an agent plays.

In these situations, the natural quantity to describe the state of agents is the distribution of agents according to their current stimuli. In large populations, the probability density function of this distribution can be approximated by a continuous function and its changes can be described by a non-linear Fokker-Planck equation. The derivation of this equation, which generalizes the replicator dynamics equation for heterogeneous populations, is contained in the Appendix of this paper. The assumptions needed for the derivation of the equation are: large population size, large number of plays of the game, and incremental (infinitesimal) structure payoffs. This approach is a well-known method for modeling multi-agent systems in the problems of physics, biology, economics and sociology, see for example, Risken (1992) and Pareschi & Toscani (2014). A similar method has been used by Traulsen, Claussen & Hauert (2005, 2006), in the analysis of the evolutionary selection by Moran process.

After deriving the equations we apply them to determine the behavior of RPS learners in symmetric 2x2 games. It should be kept in mind, however, that the time asymptotic behavior based on the Fokker-Planck equation and the time asymptotic of the original discrete-time stochastic process are not, in general, the same. Moreover, the PDE model derived in this paper is a leading approximation of a continuous-time stochastic process. Thus, any statement claimed in this paper about the convergence of the system to a particular state should be understood as a statement that the system gets close to this state within the limits of the validity of the PDE approximation.

In game theory, the method of stochastic approximation by Benaïm & Hirsch (1999) is typically used to determine the long-time behavior interacting agents. The method relies on the stability analysis of an averaged (deterministic) system of ODEs. As the dimension of a system is of the order (the number of agents ) $\times$ (the number of strategies per agent), even for small populations it is an extremely challenging task. The Fokker-Planck equation, on the other hand, is applicable for large populations, the trade-off being the loss of information about a particular player.

The analysis of the PDE model derived in this paper shows that for games with a single Nash equilibrium, the strategies of all agents converge to the dominant strategy, when the RPS rule has no memory factor, or with a memory factor and zero residuals. For learning models with a memory factor the convergence is faster than for the models with perfect memory. Additionally, the learning time in the former case varies inversely with the size of the memory factor. For games with a mixed Nash equilibrium, the learning process converges to a state in which the population mean probability equals to the equilibrium value. As the population mean probability approaches its limit, the “strength” of the learning decreases, and the individual probabilities do not, in general, converge to the equilibrium value. In this case, the population remains heterogeneous. Finally, if the memory factor is present and residuals are not zeros, agent’s strategies converge to some mixed strategy, for a generic 2x2 symmetric game.

2 The Model

We consider a series of plays of a game between randomly selected individuals in a large population. The payoff matrix of the game is given in the table 1. The behaviors are labeled A and B.

We will analyze the RPS learning in the games that a) have a single pure Nash equilibrium with $a>c,$ $d>b$ (or $a<c,$ $d<b$ ); b) two Nash equilibria with $a>c,$ $b>d$; and c) mixed Nash equilibrium when $a<c,$ $d>b.$ In the latter case, there are also two asymmetric Nash equilibria $(A,B)$ and $(B,A).$

Now, let there be a group of $N$ individuals where each individual is characterized by vector $X^{t}_{i}{}={}(S^{t}_{i,1},S^{t}_{i,2})$ representing the accumulated stimuli to play A and B, respectively, at epoch $t.$ Suppose that two agents $i$ and $j$ are selected at random to play the game. Agents play with the probabilities to cooperate $S^{t}_{1,i}/(S^{t}_{1,i}+S^{t}_{2,i})$ and $S^{t}_{1,j}/(S^{t}_{1,j}+S^{t}_{2,j}),$ for agent $i$ and agent $j,$ respectively. Denote the outcomes (A,A), (A,B), (B,A), (A,A), where the first is the action chosen by agent $i.$ As the result of the interaction agent $i$ increments his/her states according to the rule

[TABLE]

and symmetrically for agent $j.$ $(r_{1},r_{2})$ are the residual stimuli and $\mu\geq 0$ the parameter of the fading memory. The former is related to the memory factor $\bar{\mu}$ from the introduction by formula $\bar{\mu}=\mu h.$ Then, time moves to the next epoch $t+\delta$ and the process is reiterated. At time $t=0$ agents have generally different, initial levels of stimuli to cooperate and defect.

2.1 Stimuli-space

Our main interest is in the distribution of agents in the stimuli-space $x=(s_{1},\,s_{2}),$ $s_{1},\,s_{2}\geq 0,$ described by the density function (PDF) $f(x,t).$ In this space the straight lines through the origin represent the sets of stimuli of constant probability to cooperate when agent is playing a mixed strategy. The probability $m$ related to the slope $k$ as $m{}={}(1+k)^{-1}$ so that the stimuli with preference for C are located closer to the $s_{1}$-axis. For any subset $\Omega$ in the stimuli space, $\int_{\Omega}f(x,t)\,dx$ represents the proportion of agents with their stimuli in the set $\Omega$ at time $t.$

The following equation is found to be a leading order approximation for the process, when parameters $h,\delta$ are small and the number of players is large.

[TABLE]

and velocity $u$

[TABLE]

where

[TABLE]

The analysis of the equation can be understood from the behavior of the system of ODEs:

[TABLE]

and the equation for $\bar{m}(t)$ that follows from (5):

[TABLE]

The derivation of equations (5)–(8) is contained in the Appendix.

Velocity $u$ represents the rates of change of the stimuli of agents whose current state is given by $x.$ The rates are proportional to the group average payoffs for corresponding actions, and “penalized” by memory for large deviations from the default residual levels $(r_{1},r_{2}).$ Equation (8) is convenient for analyzing the games with pure Nash equilibrium. For the game with a mixed equilibrium, described by the frequency $m_{*}$ to play A, more convenient is equation

[TABLE]

The nonlinear equation (5) is the first order approximation of the stochastic, learning process. The next order approximation contains a diffusion term, with the diffusion coefficients of order $h.$ The diffusion generally prevents the convergence of learning of the group to a single strategy (fixation). For example, the convergence of the group learning to a single strategy, based on equation (5), only indicates that the distribution of agents’ strategies gets close to that particular strategy, within the limits of validity of equation (5).

Now we consider the dynamics of learning in symmetric $2\times 2$ games. First we consider models with no memory factor, $\mu=0.$

2.2 Pure Nash equilibrium

Let $a>c$ and $d>b.$ The characteristic property of this regime is the positive sign of $\bar{m}^{\prime}(t)$ in equation (8), for any distribution function $f.$ This reflects the fact that action C is your best choice, no matter what your opponent does.

It is shown in Appendix that $\bar{m}(t)$ increases to its maximum value 1 and the support of function $f$ is transported to infinity along the trajectories of ODE (7). In particular the stimuli of all agents, for large values of $t$ will be located below any straight line of positive slope through the origin. That is, asymptotically all agents learn to play the equilibrium strategy “always A”. The inclusion the second-order effects does not change the asymptotic picture.

2.3 Mixed Nash equilibrium

Let $a<c$ and $d>b.$ The equilibrium mixed strategy is $m_{*}{}={}\frac{d-b}{c-a+d-b}$ and the group average probability to play A evolves according to equation (9).

It is shown in Appendix that the dynamics of equation (5) implies that $\bar{m}(t)$ converges to the equilibrium density $m_{*}.$ The population average probability to play A asymptotically coincides with the probability $m_{*}$ at the Nash equilibrium. Unlike the pure Nash equilibrium case, in general, agents keep playing with different strategies. This can be seen from the following example. If the initial data $f_{0}(x)$ describes the population of agents playing different strategies and is such that

[TABLE]

then the trajectories of the flow generated by $u$ are straight lines and for any $t>0,$

[TABLE]

Clearly, for all $t,$ there is non-diminishing spread in the distribution function $f(x,t).$ To put it differently, there is no learning when the group average probability to play A equals $m_{*},$ because players expected payoffs from A and B are equal.

2.4 Two Nash equilibria

For this type of game $a>c$ and $d<b.$ Denote by $m_{*}{}={}\frac{b-d}{a-c+b-d}.$ As in the previous case, the group average probability to play A evolves according to equation (9). If the initial data $\bar{m}(0)>m_{*}$ velocity carries $f$ to values of stimulus values of $s_{1}$ much larger than $s_{2}.$ As the average $\bar{m}(t)$ increases, the dynamics is consistent and leads to learning of the Nash equilibrium A. With $\bar{m}(0)<m_{*}$ learning converges to the other equilibrium. In the borderline case $\bar{m}(0)=m_{*},$ is unstable: stimuli increase along the straight lines, and agents retain their initial probabilities to play A and B, but any perturbations will deviate the system to A or B equilibrium.

2.5 Memory factor RPS

Consider now models with $\mu>0.$ It can be seen from equation (6) that a sufficiently large box $[0,\hat{s}_{1}]\times[0,\hat{s}_{2}]$ is invariant under the flow of (7). This can be seen from the sign of the velocity components. We show in Appendix that with residuals $r_{1},r_{2}>0,$ RPS learning will approach an asymptotically stable point in the stimuli space, for any positive payoff rates $a,b,c,d.$ All agents will tend to play a mixed strategy $m_{*}{}={}\frac{s^{*}_{1}}{s^{*}_{1}+s^{*}_{2}}.$ When (A,A) is a Nash equilibrium, the ratio of stimuli $s^{*}_{1}/s^{*}_{2}>r_{1}/r_{2},$ and agents favor action A after learning more than at their default levels.

There is an interesting limiting case of zero residuals $r_{1}=r_{2}=0.$ For such RPS models, when the game has (A,A) as the single Nash equilibrium, all agents learn this optimal strategy, as function $f(x,t)$ converges to a delta mass supported at the point $(a/\mu,0).$

In this case, the equation for mean $\bar{m}(t)$ is closely approximated by the replicator dynamics equation (11) given below. The factor $(s_{1}(t)+s_{2}(t))^{-1}$ on the right-hand side of that equation, in this case, is of order 1. The convergence to the optimal strategy is faster than in the case of learning with perfect memory ($\mu=0$), for which the factor is of the order $(1+t)^{-1}.$

Moreover, among the models with a memory factor, the learning period appears to be shorter for larger values of $\mu,$ see figure 1 and explanations in Appendix.

Thus it appears that RPS models with small residuals and large values of $\mu$ should be preferred by natural selection. In such models agents act predominantly on the basis of the last few payoffs. In this context it is worth mentioning that one of the postulates of prospect theory of Kahneman & Trverski (1984), is the statement that people actions (in games with monetary payoffs) are directed by the increase in their total wealth, rather than the total accumulated wealth, which shows a tendency to use short memory and suggests that RPS learning might be at work.

2.6 Relation to the Replicator Dynamics Equation (RDE)

If one postulates that all agents have the same, or approximately the same, stimuli

[TABLE]

for some $s_{1}(t),s_{2}(t),$ so that $f(s,t)$ is represented by a delta function supported at $(s_{1}(t),s_{2}(t)),$ then equation (5), leads to a variant of the replicator dynamics equation for the probability to cooperate $m(t){}={}s_{1}(t)/(s_{1}(t)+s_{2}(t)):$

[TABLE]

Notice the positive factor on the right-hand side of the equation. For a learning processes in which stimuli increase the learning rate slows down. The extent to which hypothesis (10) is consistent with the dynamics of (7) is limited only to the cases when the latter has a single asymptotically stable fixed point.

3 Appendix: a PDE model

3.1 Fokker-Planck equation

Consider a group of $N$ individuals acting according to RPS learning rule described in section 2. Let $X^{t}=(X^{t}_{1},..,X^{t}_{N})$ represent the vector of pairs of stimuli for all members at epoch $t.$ Each component of this vector is 2-dimensional: $X^{t}_{i}{}={}(S^{t}_{1,i},S^{t}_{2,i}).$ By $w_{h}(\bar{x},t),$ where $x\in[0,1]^{2N},$ we denote PDF for distribution of $X^{t}.$ We will write $\bar{x}=(x_{1},..x_{N}),$ where each $x_{i}{}={}(s_{1,i},s_{2,i}).$ The probability to play A will be denoted as $\lambda_{i}{}={}s_{1,i}/(s_{1,i}+s_{2,i}).$

Suppose that member $i$ and $j$ are selected for the interaction. There will be only one game played during the period from $t$ to $t+\delta.$ The matrix of payoffs is described in table 1. The range of parameters $\delta,\,h,\,N$ will be restricted later on.

Conditioned on the event $X^{t}=\bar{x},$ the agent probabilities for the next period are set according to the RPS rule (4), which in the notation of the stochastic process are

[TABLE]

and symmetrically for $X^{t+\delta}_{j}.$ For all other agents, $X^{t+\delta}_{k}{}={}X^{t}_{k}$ for $k\not=i,j.$ The definition of $X^{t}$ makes it a discrete-time Markov process. We proceed by writing down the integral form of the Chapman-Kolmogorov equations and approximate its solution by a solution of the Fokker-Planck equation (forward Kolmogorov’s equation), for small values of $\delta,h$ and large $N.$

Change of $w_{h}(\bar{x},t)$ from $t$ to $t+\delta,$ can be described in the following way.

[TABLE]

This equation can be written in slightly different way:

[TABLE]

The above equation can be used to obtain $2N$-dimensional ODE approximation of the stochastic process by evaluating $\lim_{h\to 0}\left(\mathbb{E}[X^{t+h}\,|\,X^{t}]-X^{t}\right)/h{}={}F(X^{t}).$ This approach was implemented in the method of stochastic approximation developed by Benaim-Hirsch (1999) and applied to the study of convergence of stochastic fictitious play processes. The method guarantees the convergence of the process $X^{t}$ under certain stability conditions for the dynamics of the associated ODE.

The large dimension of that dynamical system is an obstacle for further analysis. In contrast, we would like to obtain an equation for the distribution of large number $N$ of agents in 2-dimensional stimuli space. For this, denote the PDF of the distribution by

[TABLE]

where $\bar{x}_{k}$ is a $2N-2$ dimensional vector of all coordinates, excluding $x_{k}.$ In statistical physics this function is also called one-particle distribution. In the formulas to follow we need to use two-particle distribution function

[TABLE]

where $\bar{x}_{ij}$ is the $2N-4$ dimensional vector of all coordinated excluding $x_{i}$ and $x_{j}.$ Function $g_{h}$ is symmetric in $(x,y)$ and is related to $f_{h}$ by the formulas

[TABLE]

The moments of function $f_{h}$ and $g_{h}$ are computed from the moment of $w_{h}:$

[TABLE]

and

[TABLE]

This follows from the definition of these functions.

Now we use (13) to obtain an integral equation of the change of function $f_{h}.$ For that select $\phi(\bar{x}){}={}\psi(x_{k}),$ sum over $k$ and take average. We get

[TABLE]

The right-hand side can be conveniently expressed in terms of the two-particle function $g_{h}:$

[TABLE]

where $x=(s^{x}_{1},s^{x}_{2}),$ $y=(s^{y}_{1},s^{y}_{2}),$ and $\lambda(x){}={}s^{x}_{1}/(s^{x}_{1}+s^{x}_{2}),$ and similar for $\lambda(y).$ In the processes with large number of agents and random binary interactions, two-particle distribution function can be factored into two independent distributions:

[TABLE]

With this relation, (15), becomes a family of non-linear integral relations for the next time step distribution $f_{h}(x,t+\delta).$ Taking the Taylor expansions up to the first order for the increment of the test function $\psi,$ we obtain integral equations:

[TABLE]

where $\bar{m}_{h}(t){}={}\int\lambda(x)f_{h}(x,t)\,dx.$

Combining various terms on the right-hand side of the equation we get

[TABLE]

with

[TABLE]

[TABLE]

By integrating the right-hand side by parts, and assuming that $h,\delta$ are small and $N$ is large, in such a way that $\frac{4h}{\delta N}=1,$ we obtain the Fokker-Planck equation:

[TABLE]

where $x=(s_{1},s_{2}),\,s_{1},s_{2}>0$ and the drift velocity is given by the formula

[TABLE]

where

[TABLE]

Equation (8) is obtained from (16) by multiplying it by $\lambda(x){}={}s_{1}(s_{1}+s_{2}),$ and integrating by parts. We’re assuming here that the support of $f$ is contained in the interior of the first quadrant, so that $f$ is zero on the boundary. That is to say that all agents play mixed strategies. This is a natural hypothesis, as nothing else can be learned if an agent chooses an action with certainty.

Consider now the learning from playing a symmetric $n\times n$ game. Let the payoff to playing $i^{th}$ action against $j^{th}$ be equal to $a_{ij}h.$ Denote the stimulus vector $x=(s_{1},..,s_{n}),$ and the population average probability to play $i$ by

[TABLE]

The first-order approximation of the RPS learning process is given by the Fokker-Planck equation

[TABLE]

on the domain with $s_{i}\geq 0,$ $i=1..n.$ In this equation, the velocity vector $u=(u_{1},..,u_{n})$ is given by its components

[TABLE]

Now we consider in some detail the learning in 2x2 games. Much of the analysis of equation (16) is derived from the behavior of trajectories of ODE (7). The solution $f$ of (16) is obtained by transporting the support of $f_{0}$ along trajectories of the dynamical system (7) and changing the values $f_{0}$ so that the “mass” (measured by the density function $f$) of any “fluid element” remains constant. In fact on can write down the formula for $f$ in terms of $u$ and prove that the solution $f$ of the non-linear problem exists and unique. This can be done by standard methods of PDEs, but it is outside of the scope of the present paper. Here, we will be interested in the long time, qualitative asymptotic for $f(x,t).$

Equation (16) is considered in the first quadrant $s_{1},s_{2}>0.$ For $\mu=0,$ the boundary of the domain in invariant under the flow of (7). For the model with fading memory, $\mu>0,$ the velocity $u$ at the boundary is directed into the flow domain. In either case, we will assume that the function $f_{0}(x)$ is zero on the boundary. Then, this property will hold for all $t>0.$ Additionally, in all of the analysis below, we assume that $f_{0}$ is a continuously differentiable function with compact support (zero outside some bounded set).

Consider the case of the pure Nash equilibrium ($a>c,d>b$) and no memory effects, $\mu=0.$ The velocity

[TABLE]

where matrix

[TABLE]

For any $t,$ the origin is an unstable node with two positive eigenvalues; the eigenvalue corresponding to $s_{1}$–direction is the dominant one. From the phase portrait of the ODE it is clear that the flow transports the support of $f$ into the region where $s_{1}\gg s_{2},$ which correspond to the case of all agents asymptotically in $t$ adopting choice C in the game.

In the case of the mixed Nash equilibrium ($a<c,d>b$) and no memory effect, $\mu=0,$ the origin is an unstable node. When $\bar{m}(t)=m_{*}$ then two positive eigenvalues coincide, and all trajectories are straight lines through the origin. In general, however, $\bar{m}(t)\not=m_{*},$ if for example they are not equal at time $t=0.$ In such cases equation (9) can be used to show that $\lim_{t\to\infty}\bar{m}(t)=m_{*}.$ Let $\bar{m}(0)>m_{*}.$ Then, according to (9), $\bar{m}(0)>\bar{m}(t)>m_{*}$ for all $t,$ and $\bar{m}(t)$ converges to $m_{*}$ provided that

[TABLE]

diverges. Notice also, that the derivative of ratio $s_{2}/s_{1}$ along a flow trajectory equals

[TABLE]

Thus, if at $t=0$ the support of $f_{0}(\cdot)$ is strictly inside the first quadrant, then there $c>0$ such that $s_{2}/s_{1}>c$ for all points in the support of $f(\cdot,t)$ for all later times. In particular, one can estimate

[TABLE]

Finally,since $u(x,t)$ is uniformly bounded, i.e., for any $(x,t),$ $|u(x,t)|<C,$ for some $C,$ then for any $x$ in the support of $f(\cdot,t)$ there a constant $\hat{C}$ such that $|x|<\hat{C}(1+t).$ From this and (19) it follows that

[TABLE]

for some constant $C>0,$ and so, the integral in (18) is infinite. The case $\bar{m}(0)<m_{*}$ follows by similar arguments.

Consider now the model with the memory decay when $\mu>0$ and residuals $r_{1},r_{2}>0.$ For any value of $\bar{m}(t)\in[0,1]$ and any set of positive parameters of the game $a,b,c,d>0$ the right-hand side of (7) has a steady state $(s^{0}_{1}(t),s^{0}_{2}(t))$ in the interior of the first quadrant, with $s^{0}_{1}>r_{1},$ $s^{0}_{2}>r_{2},$ and this point is an asymptotically stable node. The other steady state is the origin, which is an unstable node. The support of $f(x,t)$ moves in the direction of the stable interior point, contracting in size. When it becomes sufficiently small, the dynamics can be effectively approximated by ODE: $\frac{d(s_{1},s_{2})}{dt}{}={}(\bar{u}_{1},\,\bar{u}_{2})$ where the new velocity

[TABLE]

In a long run the fixed point will settle at the stable, interior state $x^{*}=(s^{*}_{1},s^{*}_{2}),$ and $f(x,t)$ will be a delta-function supported at that point. In this process the agents learn to play A with probability $m_{*}{}={}s^{*}_{1}/(s^{*}_{1}+s^{*}_{2}).$

A special case of zero residuals $r_{1}=r_{2}=0$ deserves a discussion. In the limit of zero residuals $r_{1},r_{2}\to 0,$ velocity $u$ (for any fixed $t>0$) has three fixed points: $x_{0}(t)=((a\bar{m}(t)+d(1-\bar{m}(t))/\mu,0),$ $x_{1}=(0,(c\bar{m}(t)+b(1-\bar{m}(t))/\mu)$ and $(0,0).$ When $a>c,d>b,$ the first, corresponding to the strategy “always A”, is an asymptotically stable node, the second, corresponding to “always B”, is a saddle, and the origin is an unstable node. One can compute that on any trajectory of the velocity field $u(x,t)$ inside the first quadrant,

[TABLE]

Thus, the population average probability to play A, $\bar{m}(t),$ increases to 1, the stable stationary point $x_{0}(t)$ converges to $x_{0}=(a/mu,0),$ and the support of $f(x,t)$ moves toward point $x_{0}.$ The agents with memory decay and zero residual levels do learn the optimal strategy. Moreover, because the learning occurs in the bounded region of the stimuli space, the convergence to the equilibrium is faster than the case of learning with perfect memory $\mu=0.$ Using equation (11) we can also estimate on the rate of convergence as a function of $\mu.$ The rate is proportional to $\mu$ implying that the characteristic time of convergence to be $C/\mu.$ Figure 1 shows the simulations of the stochastic process in this regime for different values of $\mu.$ It shows that the prediction based on the model (16) is in good agreement with the stochastic learning process. We conclude that models with low residuals and high $\mu$ perform better in learning the optimal strategy in 2x2 games and thus more likely to evolve.

3.2 Second order effects

The inclusion of the next order approximation into consideration adds a diffusion term into equation (16) with the diffusion coefficients proportional to $h.$ In the problems where the drift velocity takes $f$ to regions with large values of $x,$ as in the case of the pure or mixed Nash equilibrium, diffusion will have a marginal effect. In the problems with asymptotically stable fixed points in stimuli space (short memory models) diffusion will create a stationary distribution of $f$ near the fixed point, preventing all agents to adopt a single strategy.

The author wishes to thank the anonymous referees for patient reading of the manuscript and detailed comments that helped to improve it in so many ways.