BM model - Prisoner's Dilemma

Program to assess the impact of floating-point errors

Luis R. Izquierdo
The Macaulay Institute

In this model, players decide stochastically whether to cooperate or defect. Each player’s strategy is defined by the probability of undertaking each of the two actions available to them. After every player has selected an action according to their probabilities, every player receives the corresponding payoff and revises her strategy. The revision of strategies takes place following a reinforcement learning approach: players increase their probability of undertaking a certain action if it led to payoffs above their aspiration level, and decrease this probability otherwise. When learning, players in the BM model use only information concerning their own past choices and payoffs, and ignore all the information regarding the payoffs and choices of their counterparts. More precisely, the updating of a strategy takes place in two steps. First, each player i calculates her stimulus sai for the action just chosen a (either Cooperate (C) or Defect (D)), according to the following formula (where every variable is indexed in i):

sa = \(TraditionalForm\`\(π\_a - A\)\/sup[\(\(|\)\(T - A\)\(|\)\)\ , \ \(\(|\)\(R - A\)\(|\)\)\ , \(\(|\)\(P - A\)\(|\)\), \(\(|\)\(S - A\)\(|\)\)]\)

where (Subscript[π, a]) is the payoff obtained having selected action a, A is the player’s aspiration level, and T , R , P , S are the possible payoffs the player might receive (T for Temptation, R for Reward, P for Punishment, and S for Suckers). Hence the stimulus is always a number in the interval [-1, 1]. Secondly, having calculated their stimulus sai , each player i updates her probability pai of undertaking the selected action a as follows (where every variable is indexed in i):

(Subscript[p, a, n + 1])={{  ({{({{(Subscript[p, a,n ]) +  
                        l (Subscript[s, a, n]) (1  -  
                              (Subscript[p, a,n ])), if  (Subscript[s, a, n])}})≥0}, {({{(Subscript[p, a,n ]) +  
                        l (Subscript[s, a, n])  (Subscript[p, a,n ]),                if  (Subscript[s, a, n])}})<0}})

where a∈{C,D},  pa,n is the probability of undertaking action a in time-step n, sa,n is the stimulus experienced after having selected action a in time-step n, and l is the learning rate (0 < l < 1). Thus the higher the stimulus (or the learning rate), the larger the change in probability. The updated probability for the action not selected derives from the constraint that probabilities must add up to one.

The output of this program is a series of histograms. For each element j of the vector reportingList (see section Parameters), two histograms are created. Each of those histograms shows the probability to cooperate for one agent after j iterations calculated over runs (see section Parameters) number of runs. The difference between the two histograms is that one corresponds to simulation runs using floating-point arithmetic and the other one corresponds to simulation runs using rational arithmetic. All the histograms are saved as files in potentially various formats (e.g. gif, tiff, and eps) and different sizes.

IMPORTANT: This program assumes that both agents' probabilities change in the same way, so calculations are performed only for one of them. This assumption is only valid for very specific sets of parameters, and is made here to increase the speed of simulations. Changing the program to allow for different probabilities is straightforward.

Licence

BM model - Prisoner's Dilemma
Copyright (C) 2005 Luis R. Izquierdo

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 2
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.

Setting the Environment

Needs["Graphics`Graphics`"];
SetDirectory["C:\Temp"];

Parameters

runs = 1000;
reportingList = {1, 10, 100, 500, 1000, 
          2000};

  
\(\(SeedRandom[
        0];\)

 (*Number\ of\ categories\ in\ the\ histograms*) \)
categories  =  
        40;

 (*Initial Probability*) 
initialProb = 1/2;
\(\(initialProbFloat = 
        N[initialProb];\)

 (*Payoffs*) \)
T = 4;
R = 3;
P = 1;
\(\(S = 
        0;\)

 (*Aspiration\ Threshold*) \)
\(\(A = 
        2;\)

 (*Learning\ rate*) \)
l = 
        1/2;

  
extensions = {".gif"};
figureSize = 600;

Definition of functions and labels

\(\(\( (*Labels . \ 
          This\ is\ done\ for\ efficiency\ and\ code\ readability*) \)\(
\)\(\(Cooperate = 0;\)

    \(Defect = 1;\)

    \(CC = 0;\)

    \(CD = 1;\)

    \(DC = 2;\)

    \(DD = 
        3;\)

 (*Calculation\ of\ possible\ stimuli . \ 
          These\ operations\ are\ performed\ only\ once\ for\ efficiency*) 

    \(divisor = Max[Abs[T - A], Abs[R - A], Abs[P - A], Abs[S - A]];\)

    \(stimulusAfterCC\  = \ \((R - A)\)/divisor;\)

    \(stimulusAfterCD\  = \ \((S - A)\)/divisor;\)

    \(stimulusAfterDC\  = \ \((T - A)\)/divisor;\)

    \(stimulusAfterDD\  = \ \((P - A)\)/
          divisor;\)

 (*Stimulus\ after\ a\ certain\ outcome . \ 
          This\ function\ returns\ the\ stimulus\ after\ a\ given\ outcome\ \((e . g . \ 
            CD)\)\ for\ the\ agent\ who\ has\ selected\ the\ action\ indicated\ by\ the\ first\ letter\ of\ the\ outcome\ \(\((e . g . \ 
              C)\)\(.\)\)*) 

    \(stimulusAfter[x_] := 
        Switch[x, 
CC, 
          stimulusAfterCC, 
CD, 
          stimulusAfterCD, 
DC, 
          stimulusAfterDC, 
DD, 
          stimulusAfterDD];\)

 (*Changes\ in\ the\ probability\ of\ cooperating\ after\ given\ outcomes . \ 
          The\ first\ argument\ is\ the\ outcome, \ 
      and\ the\ second\ argument\ is\ the\ probability\ to\ be\ updated . \ 
          This\ function\ returns\ the\ change\ on\ the\ probability\ indicated\ by\ the\ second\ argument\ after\ the\ outcome\ indicated\ in\ the\ first\ argument\ \((e . g . \ 
            CD)\)\ for\ the\ agent\ who\ has\ selected\ the\ action\ indicated\ by\ the\ first\ letter\ of\ the\ outcome\ \(\((e . g . \ 
              C)\)\(.\)\)*) 

    \(changeAfterOn[x_, 
          y_] := 
If[\((x == CC)\) || \((x == 
                CD)\), \ 
 (*Agent\ has\ cooperated*) 
If[
            stimulusAfter[x] ≥ 0, \((1 - y)\)*l*
              stimulusAfter[x], 
            y*l*stimulusAfter[
                x]], 
 (*Agent\ has\ defected*) 
If[
            stimulusAfter[x] ≥ 
              0, \(-\((y*l*stimulusAfter[x])\)\), \(-\((1 - y)\)\)*l*
              stimulusAfter[x]]
];\)\)\)\)

Initialisation

noOfReports = Length[reportingList];
iterationsPerRun = 
        reportingList[[noOfReports]];
control  =  
        Ceiling[runs / 
            noOfReports];

  
decisions = {, };
\(\(decisionsFloat = {, };\)

 (*Initialisation\ of\ the\ list\ of\ histograms*) \ \)
histoList = 
        Table[0, {noOfReports}, {categories}];
histoListFloat = 
        Table[0, {noOfReports}, {categories}];

 (*List to store the average of floating - 
        point errors*) 
averageErrorSeries = 
        Table [(i - 1)*(2 - j), {i, iterationsPerRun + 1}, {j, 
            2}];
maximumErrorSeries = 
        Table [(i - 1)*(2 - j), {i, iterationsPerRun + 1}, {j, 
            2}];
finalProbSeries = 
        Table [(i - 1)*(2 - j), {i, runs}, {j, 2}];

Actual program

\(\(For[j = 1, 
        j ≤ 
          runs, \(j++\), 

If[
          Mod[j, control] == 0, Print["Starting run ", \ j]; 
          Print[Date[]];]; 

prob1 = 
          initialProb; 
prob1Float = 
          initialProbFloat; 
averageErrorSeries[\([1, 
              2]\)] = prob1 - 
            Rationalize[prob1Float, 
              0]; 
maximumErrorSeries[\([1, 2]\)] = 
          prob1 - Rationalize[prob1Float, 
              0]; 

histoCounter = 
          1; 

For[i = 1, 
          i ≤ 
            iterationsPerRun, \(i++\), 

firstRandomNumber = 
            Random[]; 
If[firstRandomNumber < prob1, 
            decisions[\([1]\)] = Cooperate, 
            decisions[\([1]\)] = Defect]; 
If[
            firstRandomNumber < prob1Float, 
            decisionsFloat[\([1]\)] = Cooperate, 
            decisionsFloat[\([1]\)] = 
              Defect]; 
 (*Remember\ that\ it\ is\ assumed\ that\ agents'\ probabilities\ are\ the\ same*) 
secondRandomNumber = 
            Random[]; 
If[secondRandomNumber < prob1, 
            decisions[\([2]\)] = Cooperate, 
            decisions[\([2]\)] = Defect]; 
If[
            secondRandomNumber < prob1Float, 
            decisionsFloat[\([2]\)] = Cooperate, 
            decisionsFloat[\([2]\)] = 
              Defect]; 

If[
            decisions\  ≠ \ decisionsFloat, 
            Print["Wrong decision taken because of floating-point numbers!"]; 
            Abort[];]; 

 (*We\ only\ calculate\ the\ new\ probability\ for\ one\ player\ because\ it\ is\ assumed\ that\ both\ probabilities\ are\ always\ the\ same*) 
outcome = 
            Switch[decisions, 
{Cooperate, Cooperate}, 
              CC, 
{Cooperate, Defect}, 
              CD, 
{Defect, Cooperate}, 
              DC, 
{Defect, Defect}, 
              DD]; 
prob1 = 
            prob1 + changeAfterOn[outcome, 
                prob1]; 

outcomeFloat = 
            Switch[decisionsFloat, 
{Cooperate, Cooperate}, 
              CC, 
{Cooperate, Defect}, 
              CD, 
{Defect, Cooperate}, 
              DC, 
{Defect, Defect}, 
              DD]; 
prob1Float = 
            prob1Float + 
              changeAfterOn[outcomeFloat, 
                prob1Float]; 

 (*Let' 
              s\ calculate\ the\ error*) 
prob1FloatDigits = RealDigits[prob1Float, 2]; 
prob1FloatRational = 
            FromDigits[prob1FloatDigits[\([1]\)], 2]/
              2^\((Length[prob1FloatDigits[\([1]\)]] - 
                    prob1FloatDigits[\([2]\)])\); 
error = 
            prob1 - prob1FloatRational; 
 (*\(Print["Error: ", 
                error];\)*) 

averageErrorSeries[\([i + 1, 
                2]\)] = \((\((j - 1)\)*averageErrorSeries[\([i + 1, 2]\)] + 
                  error)\)/j; 
If[
            Abs[error] > Abs[maximumErrorSeries[\([i + 1, 2]\)]], 
            maximumErrorSeries[\([i + 1, 2]\)]\  = \ 
              error]; 

 (*Report*) 
If[
            i == 
              reportingList[\([histoCounter]\)], 
If[
              Ceiling[prob1*categories] ≤ 
                0, \(histoList[\([histoCounter, 
                    1]\)]++\), \(histoList[\([histoCounter, 
                    Ceiling[prob1*categories]]\)]++\)]; 
If[
              Ceiling[prob1Float*categories] ≤ 
                0, \(histoListFloat[\([histoCounter, 
                    1]\)]++\), \(histoListFloat[\([histoCounter, 
                    Ceiling[
                      prob1Float*
                        categories]]\)]++\)]; 
 (*\(Print["Reporting after iteration ", \ i, \ " for histogram ", \ 
                  histoCounter];\)*) 
If[
              histoCounter < 
                noOfReports, \(histoCounter++\)];
]; 

 (*\(Print["End of iteration ", 
                i];\)*) 
]; 

finalProbSeries[\([j, 2]\)] = 
          prob1FloatRational; 
 (*\(Print["End of run ", \ j];\)*) 
];\)\)

Writing figures to file

\(\(\( (*Create\ a\ new\ Directory*) \)\(
\)\(\(newDirectory = 
        StringJoin["./l", ToString[Numerator[l]], "-", 
          ToString[Denominator[l]]];\)

    \(newDirectory = StringJoin[newDirectory, \ "_A", ToString[A]];\)

    \(newDirectory = 
        StringJoin[newDirectory, \ "_ip", 
          ToString[N[initialProb]]];\)


    \(CreateDirectory[newDirectory];\)

    \(SetDirectory[newDirectory];\)

    \(numberOfExtensions = \ \ Length[
          extensions];\)


    \(histoCategories = 
        Table[k/categories, {k, 0, categories}];\)


    \(For[i = 1, 
        i ≤ 
          noOfReports, \(i++\), 

 (*Make\ it\ proportions*) 
histoList[\([i]\)] /= 
          runs; 
histoListFloat[\([i]\)] /= 
          runs; 

Print["Histogram after iteration ", \ reportingList[\([i]\)], \ " calculated over ", \ 
          runs\ , \ " runs"]; 
fileName = \ 
          StringJoin[StringDrop[newDirectory, 2], "_i", 
            ToString[\ reportingList[\([i]\)]], "_r", 
            ToString[runs]]; 
fileNameFloat = \ 
          StringJoin[
            fileName, "_f"]; 

For[
          j = 1, \ j ≤ 
            numberOfExtensions, \(j++\), 
Export[
            StringJoin[fileName, extensions[\([j]\)]], 
            Graphics`Graphics`Histogram[histoList[\([i]\)], 
              FrequencyData -> True, 
              HistogramCategories -> histoCategories], 
            ImageSize -> figureSize]; 
Export[
            StringJoin[fileNameFloat, extensions[\([j]\)]], 
            Graphics`Graphics`Histogram[histoListFloat[\([i]\)], 
              FrequencyData -> True, 
              HistogramCategories -> histoCategories], 
            ImageSize -> 
              figureSize];
];

];\)


    \(finalProbSeries[\([All, 2]\)] = 
        1 - finalProbSeries[\([All, 2]\)];\)

    \(For[j = 1, \ 
        j ≤ numberOfExtensions, \(j++\), 
Export[
          StringJoin["finalProb", extensions[\([j]\)]], 
          ListPlot[finalProbSeries, PlotRange -> All, 
            FrameLabel -> {"Run number", "Final value"}, 
            Frame -> True, 
            TextStyle -> {"FontFamily" -> "Arial", "FontSize" -> 15}, PlotStyle -> PointSize[0.015]], 
          ImageSize -> figureSize]; 
Export[
          StringJoin["avgError", extensions[\([j]\)]], 
          ListPlot[averageErrorSeries, PlotRange -> All, 
            FrameLabel -> {"Iteration", "Average difference"}, 
            Frame -> True, 
            TextStyle -> {"FontFamily" -> "Arial", "FontSize" -> 15}, PlotStyle -> PointSize[0.015]], 
          ImageSize -> figureSize]; 
Export[
          StringJoin["maxError", extensions[\([1]\)]], 
          ListPlot[maximumErrorSeries, PlotRange -> All, 
            FrameLabel -> {"Iteration", "Maximum  difference"}, 
            Frame -> True, 
            TextStyle -> {"FontFamily" -> "Arial", "FontSize" -> 15}, PlotStyle -> PointSize[0.015]], 
          ImageSize -> 
            figureSize];
];\)


    \(Print[Date[]];\)\)\)\)