Techniques to understand computer simulations: Markov chain analysis. 
Figure 14

Luis R. Izquierdo, Segismundo S. Izquierdo, José M. Galán & José I. Santos

Program to create figures showing the transient distribution of the propensities to cooperate in a simple model where two reinforcement learners play a Prisoner's Dilemma game. The height of each patch denotes the relative probability that the system is on the state represented by the patch in any given time-step.

Licence

Copyright (C) 2008 Luis R. Izquierdo, Segismundo S. Izquierdo, José M. Galán & José I. Santos

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110 - 1301, USA.

Contact information :
Luis R. Izquierdo
University of Burgos, Spain.
e - mail : luis@izquierdo.name

Setting the Environment

SetDirectory[NotebookDirectory[]];

Parameters

tremblingHandsNoise=0;
xDimension=17;
xMin=0;
yDimension=17;
yMin=0;

(* Timesteps for which we want a figure showing the transient distribution. Initial conditions are always shown. *)
reportingList=Join[Table[i,{i,1,39}],Table[i,{i,40,100,5}],{150},{200}];

Function to create a figure from a probability distribution over the state space

numberOfDigits=Length[IntegerDigits[Last[reportingList]]];

CreateFigure[stateHistogram_List,timestep_]:=Module[{stateMatrix=N[Partition[stateHistogram,17]],my3DGraph},

AppendTo[stateMatrix,Last[stateMatrix]];
stateMatrix=Append[stateMatrix\[Transpose],Last[stateMatrix\[Transpose]]]\[Transpose];

my3DGraph=ListPlot3D[stateMatrix,PlotRange->{0,1},InterpolationOrder->0,Filling->Axis,Axes->{False,False,True},ViewPoint->{1,-1,0.5},Boxed->False,Axes->False,Mesh->None,ColorFunction->(ColorData["ThermometerColors"][((#3)^(1/10))]&),
Epilog->Inset[Framed[PaddedForm[timestep,numberOfDigits],Background->LightYellow],{1,0.7},BaseStyle->"Title",Alignment->Right]];

Export[StringJoin[ToString[timestep],".gif"],my3DGraph,ImageSize->800];
];

Creation of the transition matrix

TRANSLATION: state < - > coordinates

offset=-(xMin+xDimension*yMin)+1;
xMax=xMin+xDimension-1;
yMax=yMin+yDimension-1;
state[{x_,y_}]:=x+xDimension*y+offset;
coordinates[s_]:={Mod[s-1,xDimension]+xMin,Floor[(s-1)/xDimension]+yMin};
xCoordinate[s_]:=Mod[s-1,xDimension]+xMin;
yCoordinate[s_]:=Floor[(s-1)/xDimension]+yMin;

TRANSITIONS

(* In coordinates *)
MoveTowardsCooperation[{x_,y_}]:={Min[x+1,xMax],Min[y+1,yMax]};
MoveTowardsDefection[{x_,y_}]:={Max[x-1,xMin],Max[y-1,yMin]};

(* In states *)
CooperationState[s_]:=state[MoveTowardsCooperation[coordinates[s]]];
DefectionState[s_]:=state[MoveTowardsDefection[coordinates[s]]];

MOVEMENT RULES

MovementRules[s_]:=Module[{probCoopP1=yCoordinate[s]/(yDimension-1),probCoopP2=xCoordinate[s]/(xDimension-1),probMoveTowardsCooperation},
probMoveTowardsCooperation=((1-tremblingHandsNoise)*probCoopP1+tremblingHandsNoise*(1-probCoopP1))*((1-tremblingHandsNoise)*probCoopP2+tremblingHandsNoise*(1-probCoopP2));
Return[{Rule[{s,CooperationState[s]},probMoveTowardsCooperation],Rule[{s,DefectionState[s]},1-probMoveTowardsCooperation]}];
];

FINAL TRANSITION MATRIX

transitionMatrix=SparseArray[Flatten[Map[MovementRules,Range[xDimension*yDimension]]],{xDimension*yDimension,xDimension*yDimension}];

Transient Distributions

numberOfSteps=Last[reportingList];
reportCounter=1;

d1=StringJoin[NotebookDirectory[],"fig14"];
CreateDirectory[d1];
SetDirectory[d1];

(* initial state *)
timestep=0;
x=Table[0,{xDimension*yDimension}];
x[[state[{10,14}]]]=1;
CreateFigure[x,timestep++];

(* transient distributions *)
While[timestep≤numberOfSteps,

x=(x.transitionMatrix);

If[(timestep==reportingList[[reportCounter]]),
reportCounter++;
CreateFigure[x,timestep];
];

timestep++;
]

Asymptotic Distribution

x=x.MatrixPower[N[transitionMatrix],100000];
CreateFigure[x,∞];