# Evolutionary Games in Stochastic Environments library(MASS) ############################### # Evolutionary games function # ############################### # this function assumes there are a finite number of random variables determining # the payoff matrix and that these random variables are all log-normally distributed # INPUT # n - number of strategies # logvar.mean the log-means of the variables # logvar.cov the covariance matrix of the log variables # payoff a function of the variables that produces the payoff matrix # d fraction that die # T length of run # x0 initial condition # OUTPUT: x all frequencies over all time; rows correspend to different times EG=function(n,logvar.mean,logvar.cov,payoff,T,x0,d){ x=matrix(0,T,n) x[1,]=x0 for(t in 1:(T-1)){ var=exp(mvrnorm(1,mu=logvar.mean,Sigma=logvar.cov)) A=payoff(var) x[t+1,]=(1-d)*x[t,]+d*(diag(x[t,])%*%A%*%x[t,])/sum(x[t,]*(A%*%x[t,])) } return(x) } ############### # SAMPLE USES # ############### ################################ # Lottery model with n species # ################################ # INPUTS: # n number of competing species # fitness.min minimum log fitness of the species # fitness.max maximum log fitness of the species # sigma log-variance in fitness (same for all) # rho correlation in log-fitnesses between species # d fraction dying # T length of simulation Lottery=function(n=2,fitness.min=1,fitness.max=1.1,sigma=0,rho=1,d=0.9,T=1000){ logvar.mean=seq(fitness.min,fitness.max,length=n) logvar.cov=(1-rho)*sigma*diag(n)+rho*sigma # the covariance matrix payoff=function(var)var%o%rep(1,n) # the payoff function x0=rep(1/n,n) # initially equal frequencies of all types x=EG(n=n,logvar.mean=logvar.mean,logvar.cov=logvar.cov,payoff=payoff,T=T,x0=x0,d=d) #calling the command par(mfrow=c(1,2),mar=c(4,4.5,1,1),cex.lab=1.5,cex.axis=1.1) matplot(x,type="l",ylim=c(0,max(x)),xlab="time",ylab="frequencies") #plotting the data par(mar=c(4,0,1,1)) hist(log10(x[,2]),40,col="black",density=10,xlab=expression(paste(log[10]," frequency of least fit")),main="",yaxt="n",ylab="") } ################## # Hawk dove game # ################## # INPUTS: # V.mean mean payoff for a victory # sigma log-variance in payoff for a victory # C cost of a fight # Base base payoff for everyone # d fraction dying # T length of simulation HawkDove=function(V.mean=2,C=2,Base=3,sigma=0,d=0.1,T=5000){ logvar.cov=sigma logvar.mean=log(V.mean)-sigma/2 payoff=function(x)matrix(c((x-C)/2,x, 0,x/2),2,2,byrow=TRUE)+Base x0=c(1/2,1/2) x=EG(n=2,logvar.mean=logvar.mean,logvar.cov=logvar.cov,payoff=payoff,T=T,x0=x0,d=d) par(mfrow=c(1,2),mar=c(4,4.5,1,1),cex.lab=1.5,cex.axis=1.1) matplot(x,type="l",ylim=c(0,max(x)),xlab="time",ylab="frequencies") #plotting the data par(mar=c(4,0,1,1)) hist(log10(x[,2]),40,col="black",density=10,xlab=expression(paste(log[10]," frequency of least fit")),main="",yaxt="n",ylab="")} ###################### # Rock Paper Scissor # ###################### # INPUTS: # B.mean mean payoff for a win # B.sigma log-variance for a win # rho correlation between strategy payoffs for a win # C cost for a loss (needs to less than base payoff) # a.mean mean base payoff for everyone # a.sigma log-variance base payoff for everyone # rho correlation between log-wins and log-base # d fraction dying # T length of simulation RPS=function(B.mean=0.5,B.sigma=0,rho=0,C=1,a.mean=2,a.sigma=0,d=0.1,T=5000){ logvar.cov=matrix(0,6,6) logvar.cov[1:3,1:3]=diag(c(B.sigma,B.sigma,B.sigma)) logvar.cov[4:6,4:6]=diag(c(a.sigma,a.sigma,a.sigma)) logvar.cov[1:3,4:6]=diag(rep(rho*sqrt(B.sigma*a.sigma),3)) logvar.cov[4:6,1:3]=diag(rep(rho*sqrt(B.sigma*a.sigma),3)) logvar.mean=c(rep(log(B.mean)-B.sigma/2,3),rep(log(a.mean)-a.sigma/2,3)) payoff=function(x)matrix(c(x[4],x[4]-C,x[4]+x[1],x[5]+x[2],x[5],x[5]-C,x[6]-C,x[6]+x[3],x[6]),3,3,byrow=TRUE) x0=c(0.4,0.3,0.3) x=EG(n=3,logvar.mean=logvar.mean,logvar.cov=logvar.cov,payoff=payoff,T=T,x0=x0,d=d) par(mfrow=c(1,2),mar=c(4,4.5,1,1),cex.lab=1.5,cex.axis=1.1) matplot(x,type="l",ylim=c(0,max(x)),xlab="time",ylab="frequencies") #plotting the data par(mar=c(4,0,1,1)) hist(log10(x[,2]),40,col="black",density=10,xlab=expression(paste(log[10]," frequency of least fit")),main="",yaxt="n",ylab="") }