Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# adambear91/AA_BearRand_master_file.m

Last active Dec 22, 2020
Code for Bear Rand (2016) PNAS
 % % % % % % % % % % % % % % % % % % % % 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. % % % % % % % % % % % % % % % % % % % %%%%%%%%%%%%%%%%%%%%%%% %% DEFINE PARAMETERS %% %%%%%%%%%%%%%%%%%%%%%%% 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% STEADY STATE CALCULATION %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 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 payoffTable=zeros(M); % create payoff table for i=1:M for j=1:M payoffTable(i,j)=calculate_payoff(strats(i,:),strats(j,:),b,c,d,p); end end % run steady state using this payoff table and store result in ssDists ssDists(:,p==ps_steadyState)=steady_state_calc(payoffTable,N,w,a); % 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)... .*strats); end %%% Plot Results %%% figure plot(ps_steadyState,avgStrats); axis([0 1 0 1]) legend('sR','s1','sI','T'); xlabel('Probability of Repeated Game p') ylabel('Average Value') %%%%%%%%%%%%%%%%%%%%%%%%%%%% %% AGENT-BASED SIMULATION %% %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% 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,:) = ... agent_based_sim(N,b,c,d,p,w,mu,nGens,a); end %%% Plot Results %%% figure(2) plot(ps_agentSim,avgStrats,'--o'); axis([0 1 0 1]) legend('sR','s1','sI','T'); 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 % % INPUTS % 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) % % OUTPUTS % 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] params=[N,b,c,d,p,w,mu,nGens,a]; % Otherwise, use the initial condition you were passed else pop=popInit; params=[N,b,c,d,p,w,mu,nGens,a,popInit]; 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!') pause N=length(pop(:,1)); end end % Initialize the matrix that will store the payoffs that result from % each pair of agents interacting pairwisePayoffs=zeros(N); for i=1:N for j=i:N pairwisePayoffs(i,j)=calculate_payoff(pop(i,:),pop(j,:),b,c,d,p); pairwisePayoffs(j,i)=calculate_payoff(pop(j,:),pop(i,:),b,c,d,p); end end % (Leaves diagonal as 0s, since agents dont interact with themselves) % Initialize the average strategy vector for this simulation run avgStrat=zeros(1,4); %%%%%% % Main simulation loop for gen=1:nGens % Display every 100,000th generation in command window if mod(gen,100000)==0 gen end % 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 changeOccurred=0; % Did a mutation happen? if randrand)==1,1); value=input(index); end end
 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] % STRATEGY PARAMETERS % % 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 % ENVIRONMENT/GAME PARAMETERS % % 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 % OUTPUTS % % payoff is s1's payoff %%%%%%%%%%%%%%%%%%%%%%% %% CALCULATE PAYOFFS %% %%%%%%%%%%%%%%%%%%%%%%% % 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); end
 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 % % INPUTS % 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) % % OUTPUTS % 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 numStrats=length(payoffTable); transitionPs=zeros(numStrats); % 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); end end end % 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 transitionPs(i,i)=1-sum(transitionPs(:,i)); end % Now that you've got the transition matrix, find the eigenvector with % eigenvalue 1, and that's your steady state distribution! [vec,val]=eig(transitionPs); [~, 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! ssDist=vec(:,idxofEigWithVal1)./sum(vec(:,idxofEigWithVal1)); end 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 rhoInv=1; % 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 thisProd=1; % 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 thisProd=thisProd*(gi/fi); end % now that we've finished with the product for k, add it onto the sum rhoInv=rhoInv+thisProd; end % now that we've finished with the sum, take the inverse rhoBA=1/rhoInv; end
to join this conversation on GitHub. Already have an account? Sign in to comment