import numpy as np
from scipy import sparse
from matplotlib import pyplot as plt
from matplotlib.patches import Rectangle
%matplotlib inline
#init
def open_bar(n,x_mean):
strategy = np.zeros(n,dtype=int)
strategy[int(x_mean*n):] = 1
return strategy
def pre_game_all(n, alpha_mean=0.5):
alpha = np.ones(n)*0.5
return alpha
def pre_game_half(n, alpha_mean=0.5, x_mean=0.5, diff=0.05):
alpha = np.ones(n) - diff
alpha[int(n*x_mean*alpha_mean):int(n*x_mean)] = diff
alpha[int(n*x_mean+n*(1-x_mean)*alpha_mean):] = diff
return alpha
def pre_game_beta(n, alpha_mean=0.5, tot=1):
alpha = np.random.beta(alpha_mean*tot,(1-alpha_mean)*tot,size=n)
return alpha
def buy_round(G1, G2, strategy, alpha, kappa=1, mu=0.5):
n = len(strategy)
#calculate payoffs
n_defectors = strategy.sum()
n_cooperators = n - n_defectors
#calculate times each individual plays drunk and sober against coop and defect
#(note must remove themselves from the count)
times_drunk_coop = np.random.binomial(n_cooperators-(1-strategy),alpha)
times_drunk_def = np.random.binomial(n_defectors-strategy,alpha)
times_sober_coop = n_cooperators-(1-strategy) - times_drunk_coop
times_sober_def = n_defectors-strategy - times_drunk_def
#payoff without noise
payoff = G1[strategy,0]*times_sober_coop + G1[strategy,1]*times_sober_def
payoff += G2[strategy,0]*times_drunk_coop + G2[strategy,1]*times_drunk_def
#update alpha
alpha += kappa*alpha*(1-alpha)*(n_cooperators/n - mu)
#update strategy
partner = np.random.randint(n,size=n)
p_imit = 1/(1+np.exp(-(payoff[partner] - payoff)/n))
imitate = np.random.binomial(1,p_imit)
strategy = strategy*(1-imitate) + strategy[partner]*imitate
return strategy, alpha
def plot_drunk_prisoner(diff=0.2, B=0.5, rounds=1000, n=10000, kappa=1, alpha_mean=0.5, x_mean=0.5, mu=0.5):
G1 = np.float64([[1, -0.5],[1.5, 0]])
G2 = np.float64([[1, B],[B, 0]])
fig, axs = plt.subplots(2,2, figsize=(15,5))
((ax1,ax2),(ax3,ax4)) = axs
for c in ['k']:#,'b','c','y']:
strategy = open_bar(n,x_mean)
alpha = pre_game_half(n, alpha_mean, diff=diff)
drunk_idx = alpha>0.5
sober_idx = alpha<0.5
print(np.mean(alpha), np.mean(strategy))
x_drunk_init = (drunk_idx.sum() - strategy[drunk_idx].sum())/drunk_idx.sum()
x_sober_init = (sober_idx.sum() - strategy[sober_idx].sum())/sober_idx.sum()
for r in range(rounds):
strategy, alpha = buy_round(G2, G1, strategy, alpha, kappa=kappa)
x_drunk = (drunk_idx.sum() - strategy[drunk_idx].sum())/drunk_idx.sum()
x_sober = (sober_idx.sum() - strategy[sober_idx].sum())/sober_idx.sum()
a_drunk = np.mean(alpha[drunk_idx])
a_sober = np.mean(alpha[sober_idx])
ax1.plot(x_sober,a_sober,'b.')
ax1.plot(x_drunk,a_drunk,'r.')
x_ = (n - strategy.sum())/n
uniq_a = np.unique(alpha)
n_a = len(uniq_a)
ax3.plot(x_, np.mean(alpha),'k.')
ax2.plot(np.ones(n_a)*r,uniq_a,c+'.')
ax4.plot(r,x_drunk,'r.')
ax4.plot(r,x_sober,'b.')
ax1.plot(x_sober_init,diff,'m*',ms=14)
ax1.plot(x_drunk_init,1-diff,'m*',ms=14)
ax1.plot(x_sober,a_sober,'ys',ms=10)
ax1.plot(x_drunk,a_drunk,'ys',ms=10)
ax3.plot(x_mean,alpha_mean,'m*',ms=14)
ax3.plot(np.mean(alpha), 1-np.mean(strategy),'ys',ms=10)
for ax in [ax1,ax2]:
plt.sca(ax)
plt.ylim(0,1)
plt.xticks([])
plt.yticks(size=20)
for ax in [ax3,ax4]:
plt.sca(ax)
plt.ylim(0,1)
plt.xticks(np.arange(0,rounds+1,500), size=20)
plt.yticks(size=20)
for ax in [ax1,ax3]:
plt.sca(ax)
plt.xlim(0,1)
ax1.set_ylabel(r'$\alpha$', size = 22)
ax2.set_ylabel(r'$\alpha$', size = 22)
ax3.set_ylabel(r'$\alpha$', size = 22)
ax4.set_ylabel(r'$x$', size = 22)
ax3.set_xlabel(r'$x$', size = 22)
ax4.set_xlabel(r'Round', size = 22)
print(np.mean(alpha), np.mean(strategy))
print(a_sober,a_drunk)
plot_drunk_prisoner(diff=0.2, B=0.5, rounds=3000)
plot_drunk_prisoner(diff=0.4, B=0.5)
plot_drunk_prisoner(diff=0.5, B=0.5)
plot_drunk_prisoner(diff=0.1, B=0.5)
plot_drunk_prisoner(diff=0.5, B=0.2)
plot_drunk_prisoner(diff=0.1, B=0.2)
plot_drunk_prisoner(diff=0.5, B=0.75)
plot_drunk_prisoner(diff=0.1, B=0.75)
plot_drunk_prisoner(diff=0.4, B=0.75)