###################################################################### ##IPD: Save this file as IPD # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read IPD # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: April 2013 print(`Created: April 2013`): print(` This is IPD `): print(`To implement and illustate the Iterated Prisoner Dilemma, inspired by`): print(`the mini-revolutionary article by William H. Press`): print(`and Freeman Dyson`): print(`Iterated Prisoner's Dilemma contains strategies that dominate`): print(`any evloutionary opponent`); print(`Written by Doron Zeilberger and his Experimental Math Spring 2013`): print(`class. `): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are: Apq, PowerM, ST `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: FindGoodp, FindGoodq, `): print(`FindGoodRat, FindGoodRatS, FindVGoodp, PDv, sX, sXf1`): print(` sY, sYf1 `): print(` `): elif nops([args])=1 and op(1,[args])=Apq then print(`Apq(p,q): the matrix in Fig. 2 of the Price-Dyson paper`): print(`in PNAS `): print(` Try Apq([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5]); `): elif nops([args])=1 and op(1,[args])=FindGoodp then print(`FindGoodp(R,S,T,P,p1,p4): the vectors p such that sY(p,q,R,S,T,P)`): print(`in terms of the parameter x`): print(`does not depend on q`): print(`Try:`): print(`FindGoodp(R,S,T,P,p1,p4);`): elif nops([args])=1 and op(1,[args])=FindGoodRat then print(`FindGoodRat(R,S,T,P,p1,c): the vectors p such that `): print(`(sX(p,q,R,S,T,P)-P)-c*(sY(p,q,R,S,T,P)-P)`): print(` is the same, regardless of q, for Numeric R,S,T,P`): print(`phrased in terms of p1`): print(`followed by the feasible range of p1`): print(`Try:`): print(`FindGoodRat(3,0,5,1,p1,2);`): elif nops([args])=1 and op(1,[args])=FindGoodRatS then print(`FindGoodRatS(R,S,T,P,p1,c): the vectors `): print(` p, phrased in terms of p1, such that`): print(`(sX(p,q,R,S,T,P)-P)-c*(sY(p,q,R,S,T,P)-P)`): print(` is the same, regardless of q, for symbolic (or numeric) R,S,T,P`): print(`Try:`): print(`FindGoodRatS(R,S,T,P,p1,2);`): elif nops([args])=1 and op(1,[args])=FindVGoodp then print(`FindVGoodp(R,S,T,P,p1,c): the vectors p such that sY(p,q,R,S,T,P)=c`): print(`regardless of q`): print(`FindVGoodp(R,S,T,P,p1,c);`): elif nops([args])=1 and op(1,[args])=PDv then print(`PDv(R,S,T,P,p1,p4): the pre-computed`): print(`probability vector in Eq. 8 of of`): print(`the Press-Dyson PNAS paper`): print(`Try:`): print(`PDv(R,S,T,P,p1,p4);`): elif nops([args])=1 and op(1,[args])=PowerM then print(`PowerM(A,m): the m-th power of the matrix A`): print(`Try:`): print(`PowerM([[1,2],[2,1]],3);`): elif nops([args])=1 and op(1,[args])=sX then print(`sX(p,q,R,S,T,P): The pay-off for player X in the Iterated`): print(`Prisoner's dilemma, Try`): print(`sX([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1);`): elif nops([args])=1 and op(1,[args])=sXf1 then print(`sXf1(p,q,R,S,T,P,m,p0,q0): The pay-off for player X in the Iterated`): print(`Prisoner's dilemma with strategies p, q, at the m-th round`): print(`if player X (Y) cooporates at the first round with probability`): print(`p0 (q0). Try:`): print(`sXf1([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1,10,1,1);`): elif nops([args])=1 and op(1,[args])=sY then print(`sY(p,q,R,S,T,P): The pay-off for player Y in the Iterated`): print(`Prisoner's dilemma, Try`): print(`sY([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1);`): elif nops([args])=1 and op(1,[args])=sYf1 then print(`sYf1(p,q,R,S,T,P,m,p0,q0): The pay-off for player Y in the Iterated`): print(`Prisoner's dilemma with strategies p, q, at the m-th round`): print(`if player X (Y) cooporates at the first round with probability`): print(`p0 (q0). Try:`): print(`sYf1([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1,10,1,1);`): elif nops([args])=1 and op(1,[args])=ST then print(`ST(A): given a stochastic matrix A finds the steady-state vector`): print(`Try:`): print(`ST([[1/3,2/3],[1/4,3/4]]);`): print(``): elif nops([args])=1 and op(1,[args])=Xoptions then print(`Xoptions(R,S,T,P,m): all the options of X to fix sY`): print(`regardless of Y's action with resolution m,`): print(`followed by Y's payoff with that option. `): print(`followed by X's smallest and largest scores followed`): print(`by Y's strategy to achieve it. Try:`): print(`Xoptions(3,0,5,1,10);`): else print(`There is no ezra for`,args): fi: end: #ST(A): given a stochastic matrix A finds the steady-state vector #Try #ST([[1/3,2/3],[1/4,3/4]]); ST:=proc(A) local i,v,eq,var,n,L,tot,j: n:=nops(A): var:={seq(v[i],i=1..n)}: eq:={seq(add(A[i][j]*v[i],i=1..n)-v[j],j=1..n)}: var:=solve(eq,var): L:=[seq(subs(var,v[i]),i=1..n)]: tot:=convert(L,`+`): [seq(normal(L[i]/tot),i=1..n)]: end: #Apq(p,q): the matrix in Fig. 2 of the Price-Dyson paper #in PNAS #Try Apq([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5]); Apq:=proc(p,q): [ [p[1]*q[1],p[1]*(1-q[1]),(1-p[1])*q[1],(1-p[1])*(1-q[1])], [p[2]*q[3],p[2]*(1-q[3]),(1-p[2])*q[3],(1-p[2])*(1-q[3])], [p[3]*q[2],p[3]*(1-q[2]),(1-p[3])*q[2],(1-p[3])*(1-q[2])], [p[4]*q[4],p[4]*(1-q[4]),(1-p[4])*q[4],(1-p[4])*(1-q[4])] ]: end: #sX(p,q,R,S,T,P): The pay-off for player X in the Iterated #Prisoner's dilemma, Try #sX([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1); sX:=proc(p,q,R,S,T,P) local v: v:=ST(Apq(p,q)): normal(R*v[1]+S*v[2]+T*v[3]+P*v[4]): end: #sY(p,q,R,S,T,P): The pay-off for player Y in the Iterated #Prisoner's dilemma, Try #sY([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1); sY:=proc(p,q,R,S,T,P) local v: v:=ST(Apq(p,q)): normal(R*v[1]+T*v[2]+S*v[3]+P*v[4]): end: #PDv(R,S,T,P,p1,p4): the probability vector in Eq. 8 of of #the Press-Dyson PNAS paper #Try: #PDv(R,S,T,P,p1,p4); PDv:=proc(R,S,T,P,p1,p4) local p2,p3: p2:=(p1*(T-P)-(1+p4)*(T-R))/(R-P): p3:=((1-p1)*(P-S)+p4*(R-S))/(R-P): normal([p1,p2,p3,p4]): end: #FindGoodp(R,S,T,P,p1,p4): the vectors p such that sY(p,q,R,S,T,P) #in terms of the parameter x #does not depend on q #Try: #FindGoodp(R,S,T,P,p1,p4); FindGoodp:=proc(R,S,T,P,p1,p4) local eq,p,i,var,fd,ku,tov,ma,q,var1,g: fd:=sY(p,[1/4$4],R,S,T,P): eq:={}: for i from 2 to 4 do g:=sY(p,[0$i,1,0$(4-i)],R,S,T,P): ku:=numer(g-fd): eq:=eq union {ku}: od: var:={solve(eq,{p[2],p[3]})}: tov:={}: for var1 in var do ma:=subs(var1,[seq(p[i],i=1..4)]): if {seq(normal(diff(sY(ma,q,R,S,T,P),q[i])),i=1..4)}={0} then tov:=tov union {ma}: fi: od: if nops(tov)>1 then print(`There are more than one solution`): RETURN(tov): fi: subs({p[1]=p1,p[4]=p4},tov[1]): end: #Xoptions(R,S,T,P,m): all the options of X to fix sY #regardless of Y's action with resolution m, #followed by Y's payoff with that option. Try: #Xoptions(3,0,5,1,10); Xoptions:=proc(R,S,T,P,m) local i,j,p1,p4,lu,gu,lu1,q,f,f1,f2,i1,g: lu:=FindGoodp(R,S,T,P,p1,p4): gu:={}: for i from 0 to m do for j from 0 to m do lu1:=subs({p1=i/m,p4=j/m},lu): if lu1<>[1,1,0,0] then if lu1[2]>=0 and lu1[2]<=1 and lu1[3]>=0 and lu1[3]<=1 then g:=sY(lu1,q,R,S,T,P): f:=sX(lu1,q,R,S,T,P): f1:=Optimization[Minimize](f,seq(q[i1]=0..1,i1=1..4)): f2:=Optimization[Maximize](f,seq(q[i1]=0..1,i1=1..4)): gu:=gu union {[lu1, g,f1,f2]}: fi: fi: od: od: gu: end: #FindGoodq(R,S,T,P,q1,q4): the vectors q such that sX(p,q,R,S,T,P) #in terms of the parameters #does not depend on p #Try: #FindGoodq(R,S,T,P,q1,q4); FindGoodq:=proc(R,S,T,P,q1,q4) local eq,p,i,var,fd,ku,tov,ma,q,var1,g: fd:=sX([1/4$4],q,R,S,T,P): eq:={}: for i from 2 to 4 do g:=sX([0$i,1,0$(4-i)],q,R,S,T,P): ku:=numer(g-fd): eq:=eq union {ku}: od: var:={solve(eq,{q[2],q[3]})}: tov:={}: for var1 in var do ma:=subs(var1,[seq(q[i],i=1..4)]): if {seq(normal(diff(sX(p,ma,R,S,T,P),p[i])),i=1..4)}={0} then tov:=tov union {ma}: fi: od: if nops(tov)>1 then print(`There are more than one solution`): RETURN(tov): fi: subs({q[1]=q1,q[4]=q4},tov[1]): end: #FindGoodRatS(R,S,T,P,p1,c): the vectors p such that #phrased in terms of p1, such that #(sX(p,q,R,S,T,P)-P)=c*(sY(p,q,R,S,T,P)-P) regardless of q #For Symbolic R,S,T,P #Try: #FindGoodRatS(R,S,T,P,p1,2); FindGoodRatS:=proc(R,S,T,P,p1,c) local eq,p,i,var,fd,tov,ma,q,var1,g,p4: fd:=(sX(p,[1/4$4],R,S,T,P)-P)-c*(sY(p,[1/4$4],R,S,T,P)-P): eq:={numer(fd)}: for i from 1 to 4 do g:=(sX(p,[0$(i-1),1,0$(4-i)],R,S,T,P)-P)- c*(sY(p,[0$(i-1),1,0$(4-i)],R,S,T,P)-P): eq:=eq union {numer(g)}: od: var:={solve(eq,{p[1],p[2],p[3],p[4]})}: tov:={}: for var1 in var do ma:=subs(var1,[seq(p[i],i=1..4)]): if normal((sX(ma,q,R,S,T,P)-P)-c*(sY(ma,q,R,S,T,P)-P))=0 then tov:=tov union {ma}: fi: od: if nops(tov)>1 then print(`There are more than one solution`): RETURN(tov): elif nops(tov)=0 then print(`There are no solution`): RETURN(tov): else RETURN(subs({p[1]=p1,p[4]=p4},tov[1])): fi: end: #FindVGoodp(R,S,T,P,p1,c): the vectors p such that sY(p,q,R,S,T,P)=c #regardless of q #FindVGoodp(R,S,T,P,p1,c); FindVGoodp:=proc(R,S,T,P,p1,c) local eq,p,i,var,fd,ku,tov,ma,q, var1,g,freeman: fd:=sY(p,[1/4$4],R,S,T,P): eq:={fd-c}: for i from 1 to 4 do g:=sY(p,[0$i,1,0$(4-i)],R,S,T,P): ku:=numer(g-c): eq:=eq union {ku}: od: var:={solve(eq,{p[2],p[3],p[4]})}: tov:={}: for var1 in var do ma:=subs(var1,[seq(p[i],i=1..4)]): if normal(sY(ma,q,R,S,T,P)-c)=0 then tov:=tov union {ma}: fi: od: if nops(tov)>1 then print(`There are more than one solution`): RETURN(tov): elif nops(tov)=0 then print(`There are no solutions`): RETURN(tov): else tov:=subs({p[1]=p1},tov[1]): freeman:=solve({seq(tov[i]>=0,i=1..4),seq(tov[i]<=1,i=1..4)},{p1}): RETURN(tov,freeman): fi: end: #FindGoodRat(R,S,T,P,p1,c): the vectors p such that #(sX(p,q,R,S,T,P)-P)=c*(sY(p,q,R,S,T,P)-P) regardless of q #For numeric R,S,T,P #Try: #FindGoodRat(R,S,T,P,p1,2); FindGoodRat:=proc(R,S,T,P,p1,c) local eq,p,i,var,fd,tov,ma,q, var1,g,gu,lu,p4: fd:=(sX(p,[1/4$4],R,S,T,P)-P)-c*(sY(p,[1/4$4],R,S,T,P)-P): eq:={numer(fd)}: for i from 1 to 4 do g:=(sX(p,[0$(i-1),1,0$(4-i)],R,S,T,P)-P)- c*(sY(p,[0$(i-1),1,0$(4-i)],R,S,T,P)-P): eq:=eq union {numer(g)}: od: var:={solve(eq,{p[1],p[2],p[3],p[4]})}: tov:={}: for var1 in var do ma:=subs(var1,[seq(p[i],i=1..4)]): if normal((sX(ma,q,R,S,T,P)-P)-c*(sY(ma,q,R,S,T,P)-P))=0 then tov:=tov union {ma}: fi: od: if nops(tov)>1 then print(`There are more than one solution`): RETURN(tov): elif nops(tov)=0 then print(`There are no solution`): RETURN(tov): else gu:=subs({p[1]=p1,p[4]=p4},tov[1]): lu:=solve( {seq (gu[i]>=0, i=1..nops(gu)),seq (gu[i]<=10, i=1..nops(gu))}, {p1}): RETURN(gu,lu): fi: end: #PowerM(A,m): the m-th power of the matrix A #Try: #PowerM([[1,2],[2,1]],3); PowerM:=proc(A,m) local i,j,n,A1,k: option remember: n:=nops(A): if m=0 then RETURN([seq([0$(i-1),1,0$(n-i)],i=1..n)]): else A1:=PowerM(A,m-1): RETURN([seq([seq(add(A1[i][k]*A[k][j],k=1..n),j=1..n)],i=1..n)]): fi: end: #sXf1(p,q,R,S,T,P,m,p0,q0): The pay-off for player X in the Iterated #Prisoner's dilemma with strategies p, q, at the m-th round #if player X (Y) cooporates at the first round with probability #p0 (q0). Try: #sXf1([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1,10,1,1); sXf1:=proc(p,q,R,S,T,P,m,p0,q0) local v,A,n,i,j: v:=[p0*q0,p0*(1-q0),(1-p0)*q0,(1-p0)*(1-q0)]: A:=PowerM(Apq(p,q),m): n:=nops(A): v:=[seq(add(A[i][j]*v[i],i=1..n),j=1..n)]: R*v[1]+S*v[2]+T*v[3]+P*v[4]: end: #sYf1(p,q,R,S,T,P,m,p0,q0): The pay-off for player Y in the Iterated #Prisoner's dilemma with strategies p, q, at the m-th round #if player X (Y) cooporates at the first round with probability #p0 (q0). Try: #sYf1([1/8,1/8,1/8,5/8],[1/5,1/5,1/5,2/5],3,0,5,1,10,1,1); sYf1:=proc(p,q,R,S,T,P,m,p0,q0) local v,A,n,i,j: v:=[p0*q0,p0*(1-q0),(1-p0)*q0,(1-p0)*(1-q0)]: A:=PowerM(Apq(p,q),m): n:=nops(A): v:=[seq(add(A[i][j]*v[i],i=1..n),j=1..n)]: R*v[1]+T*v[2]+S*v[3]+P*v[4]: end: