Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Code for Bear Rand (2016) PNAS
% % % % % % % % % % % % % % % % % % %
% This MATLAB script runs the low mutation limit (steady state) calculations
% shown in Figures 2a and S1-S3, as well as the agent-based simulations
% for higher mutation rates shown in Figure S3, and plots the results.
% This script calls the following functions:
% calculate_payoff(strategy1,strategy2,b,c,d,p)
% - returns the payoff of strategy1 against strategy2 given
% the specified values of b, c, d, and p
% steady_state_calc(payoffTable,N,w,a)
% - returns the steady state distribution of the population in
% the limit of low mutation, where payoffTable is an MxM matrix
% with payoffTable(i,j) being the payoff of strategy i against
% strategy j; N is the population size; w is the selection
% strength; and a is the amount of assortment
% agent_based_sim(N,b,c,d,p,w,mu,nGens,a)
% - returns the average value of each strategy parameter for
% an agent based simulation using the specified parameter
% values, run with a mutation rate of mu and last nGens
% generations.
% % % % % % % % % % % % % % % % % % %
clear all
% Parameters used from main text, Figure 2
N = 50; b = 4; c = 1; d = 1; w = 6; a = 0;
% Agent sim parameters from Figure S3a in Supplementary Information
mu = .05; nGens = 1e7; % -> 10mil generations
% We loop through values of p (repeated game probability)
ps_steadyState = .01:.01:.99;
% Because the agent based simulations take much longer, we use a coarser
% grain
ps_agentSim = .1:.1:.9;
%%% Generate Strategy Space %%%
% For the steady state calculation, we must generate a discrete strategy
% set. To do so:
% sR, s1, and sI are all restricted to being either 0 or 1;
% we discretize the deliberation threshold (T) in intervals of d/10
% generate all possible combinations of the first three strat parameters (sR,s1,sI)
% set to be 0 or 1
[x1,x2,x3] = ind2sub([2,2,2], find(zeros(2,2,2)==0));
strats = [x1 x2 x3] - 1;
% we now combine the above with the 11 possible deliberation thresholds,
% running from 0 to d in increments of d/10
Tvals = 0:d/10:d;
strats = repmat(strats,length(Tvals),1);
Tvals = sort(repmat(Tvals,1,length(strats)/length(Tvals)));
strats = [strats Tvals'];
% ... to create a total of 2*2*2*11 (88) possible strategies
M = size(strats,1); % define M to be number of strats
%%% Main Steady State Loop %%%
% initialize matrix to store steady state distributions for each p value
ssDists = zeros(M,length(ps_steadyState));
% initialize matrix to store average strat values for each p value
avgStrats = zeros(length(ps_steadyState),4);
for p=ps_steadyState
p % display p value that's running
% initialize payoff table
% create payoff table
for i=1:M
for j=1:M
% run steady state using this payoff table and store result in ssDists
% use ssDists to get average strategy values by taking a weighted average of
% the strat values based on the distribution in ssDists
avgStrats(p==ps_steadyState,:) = sum(repmat(ssDists(:,p==ps_steadyState),1,4)...
%%% Plot Results %%%
plot(ps_steadyState,avgStrats); axis([0 1 0 1])
xlabel('Probability of Repeated Game p')
ylabel('Average Value')
%%% Main Loop %%%
% initialize matrix to store average strat values for each p value
avgStrats = zeros(length(ps_agentSim),4);
for p=ps_agentSim
p % display p value that's running
% run agent sim for select p value
avgStrats(p==ps_agentSim,:) = ...
%%% Plot Results %%%
plot(ps_agentSim,avgStrats,'--o'); axis([0 1 0 1])
xlabel('Probability of Repeated Game p')
ylabel('Average Value')
function [avgStrat,params]=agent_based_sim(N,b,c,d,p,w,mu,nGens,a,popInit)
% Agent based evolutionary dynamics simulations using the Moran process with
% exponential fitness
% N = population size
% w = selection strength
% mu = mutation rate
% nGens = number of generations per simulation
% popInit = Nx4 vector of how you want the population initialized (if this
% is left off, random initial conditions are generated)
% b,c,d,p are the game parameters used in the payoff function
% a = assortment (probability of interacting with a player having the same
% strategy as you, rather than a randomly selected player)
% avgStrat = 1x4 vector with cumulative average value of the 4 strategy parameters
% [sR,s1,sI,T]
% params = matrix of parameters used (optional)v
% Initialize population - Nx4 matrix, where
% strat consists of [sR,s1,sI,T],
% with sR,s1,sI ranging from 0 to 1 and T ranging from 0 to d
% (see calculate_payoff.m for details)
% If fewer than 9 inputs were passed to the fn (ie popInit was left
% off) start from a random initial population
if nargin<10
pop=rand(N,4); % pick random values in [0,1] for each param
pop(:,4)=pop(:,4)*d; % multiply T by d so it's in [0,d]
% Otherwise, use the initial condition you were passed
if N~=length(pop(:,1)) % just in case N and the size of popInit are different, let the user know
% and then use the N implied by the initial condition
disp('N and size of initial pop vector dont match!')
% Initialize the matrix that will store the payoffs that result from
% each pair of agents interacting
for i=1:N
for j=i:N
% (Leaves diagonal as 0s, since agents dont interact with themselves)
% Initialize the average strategy vector for this simulation run
% Main simulation loop
for gen=1:nGens
% Display every 100,000th generation in command window
if mod(gen,100000)==0
% To make the code run faster, we only want to update
% pairwisePayoffs in a given generation if the population make-up
% changed. So this variable will keep track of whether, in the
% curent generation, a change occurred. We initialize it as 0 (ie
% change did not occur) and then may change it to 1 below if a
% change is made to pop
% Did a mutation happen?
if rand<mu
% Mutation occurs:
% Pick a learner
% Pick a random strategy for her
pop(learner,:)=rand(1,4); % pick random values in [0,1] for each param
pop(learner,4)=pop(learner,4)*d; % multiply T by d so it's in [0,d]
% Indicate that a change occurred to the population
% No mutation:
% Pick a random learner
% Pick someone to imitate, proportional to fitness
% (which is exp(w*(a * payoff against own strategy + (1-a)*(avg payoff against all other strategies))))
a*diag(pairwisePayoffs)... % a * (diagonal of pairwisePayoffs)
+(1-a)*sum(pairwisePayoffs.*~eye(N),2)/(N-1)))); %(1-a) * (row-sum of all non-diagonal elements of pairwisePayoff, divided by (N-1) to make it an average)
% Check if a change will happen
% Update learner's strategy
% Now that we are done with potential mutation/learning, we want to
% update pairwisePayoffs, but only if pop may have changed (to save
% computation time)
if changeOccurred==1
% go thru each agent, and update the payoff that agent gets
% against learner, as well as payoff learner gets against that
% agent
for i=1:N
% add this generation to avgStrat
% Now we are done with the main simulation loop - this simulation run
% is over.
% convert avgStrat to an average from sum of averages across all gens
function [index,value]=pickOneProportionally(input)
% 'input' should be an Nx1 vector of values.
% this fn then picks an entry proportional to value, and
% returns that entry's index and value.
% % set NaNs to 0 - NaNs should never occur, but just
% badIdx=find(isnan(input));
% input([badIdx badIdx+length(input)/2])=0;
% first, check if input is all 0s, in which case you pick one at random
if sum(input)==0
% otherwise:
% normalize values so they sum to 1
% make a new vector 'scaled' is the CDF of input
% (i.e. that effectively lines up the values along the number line from 0 to 1)
% for example, say that you had
% input=[4 2 1 3 10]
% then after normalization you'd have input=[.2 .1 .05 .15 .5]
% and the CDF would be [.2 .3 .35 .5 1]
for i=2:length(inputN)
% now generate a random number between 0 and 1, and see which bin of scaled
% it lands in: that is, find the first entry in scaled which is larger than
% the random number
function payoff = calculate_payoff(s1,s2,b,c,d,p)
%% SETUP %%
% s1 and s2 are agent 1's and agent 2's respective strategy profiles
% where a strategy profile consists of [pR,p1,pH,T]
% pR is probability of cooperating when agent deliberates and realizes the
% game is repeated/reciprocal
% p1 is probability of cooperating when agent deliberates and realizes the
% game is 1-shot/social dilemma
% pI is probability of cooperating when agent uses intuition (and therefore
% can't discriminate repeated vs. 1-shot games)
% T is (threshold) max cost in [0,d] agent is willing to pay to deliberate
% b is benefit of cooperation
% c is cost of cooperation
% d is max possible cost of deliberation (such that the cost of
% deliberation varies uniformly in the interval [0,d])
% p is probability of repeated game
% payoff is s1's payoff
% Rename each agent's threshold to T1 and T2, respectively
T1 = s1(4);
T2 = s2(4);
% Calculate repeated game payoffs
% Repeated Game Payoff Matrix
repPayoffMat=[0 0
0 (b-c)];
% Rows correspond to agent 1's strategy
% where row 1 = agent 1 defects, row 2 = agent 1 cooperates
% Columns correspond to agent 2's strategy
% where col 1 = agent 2 defects, col 2 = agent 2 cooperates
% Values correspond to agent 1's payoffs
% e.g., if agent 1 cooperates and agent 2 defects,
% agent 1 gets a payoff of 0 because repPayoffMat(2,1) = 0
% We consider repeated-game payoffs in 4 situations:
% rPayDD: both agents deliberate
% rPayDI: only agent 1 deliberates
% rPayID: only agent 2 deliberates
% rPayII: neither agent deliberates
% E.g., to calculate rPayDI,
% We know agent 1 deliberates and realizes they're in a repeated game =>
% their prob. of cooperating is pR, i.e., s1(1)
% We know agent 1 uses their intuitive response =>
% their prob. of cooperating is pI, i.e., s2(3)
% Thus, rPayDI is a weighted average of the payoff agent 1 gets when both
% cooperate (repPayoffMat(2,2)), when only player 1 cooperates
% (repPayoffMat(2,1)), and so on
% where the prob. both cooperate is given by s1(1)*s2(3), prob. only
% player 1 cooperates is given by s1(1)*(1-s2(3)), and so on
rPayDD=s1(1)*s2(1)*repPayoffMat(2,2) + s1(1)*(1-s2(1))*repPayoffMat(2,1) + ...
(1-s1(1))*s2(1)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(1))*repPayoffMat(1,1);
rPayDI=s1(1)*s2(3)*repPayoffMat(2,2) + s1(1)*(1-s2(3))*repPayoffMat(2,1) + ...
(1-s1(1))*s2(3)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(3))*repPayoffMat(1,1);
rPayID=s1(3)*s2(1)*repPayoffMat(2,2) + s1(3)*(1-s2(1))*repPayoffMat(2,1) + ...
(1-s1(3))*s2(1)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(1))*repPayoffMat(1,1);
rPayII=s1(3)*s2(3)*repPayoffMat(2,2) + s1(3)*(1-s2(3))*repPayoffMat(2,1) + ...
(1-s1(3))*s2(3)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(3))*repPayoffMat(1,1);
% Calculate 1-shot game payoffs
% (these work the same way as the repeated game payoffs)
% 1-shot Game Payoff Matrix
onesPayoffMat=[0 b
-c (b-c)];
oPayDD=s1(2)*s2(2)*onesPayoffMat(2,2) + s1(2)*(1-s2(2))*onesPayoffMat(2,1) + ...
(1-s1(2))*s2(2)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(2))*onesPayoffMat(1,1);
oPayDI=s1(2)*s2(3)*onesPayoffMat(2,2) + s1(2)*(1-s2(3))*onesPayoffMat(2,1) + ...
(1-s1(2))*s2(3)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(3))*onesPayoffMat(1,1);
oPayID=s1(3)*s2(2)*onesPayoffMat(2,2) + s1(3)*(1-s2(2))*onesPayoffMat(2,1) + ...
(1-s1(3))*s2(2)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(2))*onesPayoffMat(1,1);
oPayII=s1(3)*s2(3)*onesPayoffMat(2,2) + s1(3)*(1-s2(3))*onesPayoffMat(2,1) + ...
(1-s1(3))*s2(3)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(3))*onesPayoffMat(1,1);
% Calculate agent 1's final payoff based on
% Prob. of repeated game (p) vs. 1-shot game (1-p)
% Prob. agents 1 and 2 deliberate (T1/d and T2/d) vs. use intuition (1-T1/d and 1-T2/d)
% Average costs agent 1 and 2 pay when deliberating (T1/2 and T2/2)
% (because the cost of deliberation varies uniformly, an agent who's
% willing to pay max cost of deliberation T pays an average cost of
% T/2 to deliberate)
% We take a weighted avg. of agent 1's payoff in 4 situations:
% 1. payoff when both agents deliberate = (p*rPayDD + (1-p)*oPayDD - T1/2)
% this happens with probability (T1/d)*(T2/d)
% 2. payoff when only agent 1 deliberates = (p*rPayDI + (1-p)*oPayDI - T1/2)
% this happens with probability (T1/d)*(1 - T2/d)
% ... and so on
payoff = (T1/d)*(T2/d)*(p*rPayDD + (1-p)*oPayDD - T1/2) + ...
(T1/d)*(1 - T2/d)*(p*rPayDI + (1-p)*oPayDI - T1/2) + ...
(1 - T1/d)*(T2/d)*(p*rPayID + (1-p)*oPayID) + ...
(1 - T1/d)*(1 - T2/d)*(p*rPayII + (1-p)*oPayII);
function [ssDist,transitionPs]=steady_state_calc(payoffTable,N,w,a)
% Calculate the steady state distribution of the Moran process with
% exponential fitness, assuming the limit of low mutation
% payoffTable is an MxM matrix of payoffs, where M is the number of strategies
% N is the population size
% w is the selection strength (0=neutral drift, ->Inf=deterministic)
% a is assortment (probability of interacting with a player having the same
% strategy as you, rather than a randomly selected player)
% ssDist is an Mx1 vector with the steady state frequency of each strategy
% transitionPs is an MxM vector with transition probilities between
% homogeneous states of the population
% To do the steady state calc, we need to assemble a transition matrix
% betewen homogeneous states of the population
for s1=1:numStrats
for s2=1:numStrats
if s1~=s2
% for each combination of strategies s1 and s2, the transition prob
% going from s1 to s2 is the product of the probability a s1 mutant
% arises, and the probability that mutant goes to fixation
%prob s1 mutant arises %prob that s1 goes to fixation
%(assuming uniform mutation)
transitionPs(s1,s2)=(1/(numStrats-1)) *calc_invasionprob(payoffTable,N,w,s1,s2,a);
% The diagonal entries of the transition matrix are the probability you
% stay in the given state - which is just 1 minus all the transition
% probabilities out of that state
for i=1:numStrats
% Now that you've got the transition matrix, find the eigenvector with
% eigenvalue 1, and that's your steady state distribution!
[~, idxofEigWithVal1]=min(abs(diag(val)-1)); % For some reason, sometimes Matlab gets the eigenvalue wrong s/t its not literally 1, just very close
% normalize the eigenvector corresponding to eigenvalue so that it sums to 1, and you're done!
function rhoBA=calc_invasionprob(payoffTable,N,w,sA,sB,a)
% here we calculate rhoAB, the probability that a sA mutant goes to fixation in
% a sB population
% this probability is the inverse of (1 + a sum of products) so I
% start by intitalizing rhoInv to 1, and then add on the sum of series of
% products, and then take the inverse
% loop through each possible number of steps the sA mutant would need to go through in order to take over
% (this ranges from 1, if there are already N-1 sA players,
% to N-1, if there is only the 1 initial sA mutant)
for k=1:N-1
% same thing - initialize the eventual product of various terms as 1,
% then multiply it by stuff
% loop through each possible make-up of the population (indexed by k=the
% number of sA players) that would need to be traversed to get from
% there being only 1 sA player to there being k sA players
for i=1:k
% for each population make-up, multiply by the ratio of the fitness
% of an sA player in the current pop to the fitness of an sB player
% in the current pop (where there are i sA players)
% fitness = exp( w * expectedpayoff)
% expectedpayoff = a * payoff against your own strategy
% + (1-a) *( (payoff against sB * fraction of other players who are sB)
% +(payoff against sA * fraction of other players who are sA)
% )
fi=exp(w* (a*payoffTable(sA,sA)+(1-a)*(payoffTable(sA,sB)*(N-i)/(N-1) + payoffTable(sA,sA)*(i-1)/(N-1))));
gi=exp(w* (a*payoffTable(sB,sB)+(1-a)*(payoffTable(sB,sB)*(N-i-1)/(N-1) + payoffTable(sB,sA)*(i)/(N-1) )));
% multiple the current state of the product by the ratio for this value of i
% now that we've finished with the product for k, add it onto the sum
% now that we've finished with the sum, take the inverse
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.