Game Theory: Parrondos Paradox

We import the libraries

In [1]:
import numpy as np
import matplotlib.pylab as plt

We define the first game

In [2]:
def GameA(Money2):
    Money=Money2.copy()
    eps=0.005
    pg=0.5-eps
    rga=np.random.random()
    if(rga<pg): 
        Money=np.append(Money,Money[-1]+1)
    else:
        Money=np.append(Money,Money[-1]-1)
    return Money

We define the second game

In [3]:
def GameB(Money2):
    Money=Money2.copy()
    eps=0.005
    p2=0.9+eps
    p3=0.25+eps
    M=Money[-1]
    ra=np.random.random()
    if (M%3==0):
        if (ra<p2):
            Money=np.append(Money,M-1)
        else:
            Money=np.append(Money,M+1)
    else:   
        if(ra<p3): 
            Money=np.append(Money,M-1)
        else:
            Money=np.append(Money,M+1)
    return Money

Let us explore how the behavior of the second game is

In [4]:
moneyB=np.array([1001])
for i in range(1000):
    r=np.random.random()
    moneyB=GameB(moneyB)
print(moneyB)

[1001 1002 1001 ...  999  998  997]
In [5]:
plt.plot(moneyB,label='B')
Out[5]:
[<matplotlib.lines.Line2D at 0x16f385c6d30>]

As it is not clear, let us see the averages, so we will have 3 arrays, one for game A, another for game 2, and the third one for the combination

In [6]:
averaging=5000
N=200
aveA=np.zeros(N)
aveB=np.zeros(N)
ave=np.zeros(N)
for j in range(averaging):
    moneyA=np.array([0])
    moneyB=np.array([0])
    money=np.array([0])
    for i in range(N-1):
        r=np.random.random()
        moneyA=GameA(moneyA)
        moneyB=GameB(moneyB)
        
        if r<0.5:
            money=GameA(money)
        else:
            money=GameB(money)
    aveA+=moneyA
    aveB+=moneyB
    ave+=money
aveA=aveA/averaging
aveB=aveB/averaging
ave=ave/averaging

And the plot of the averages is,

In [7]:
plt.xlabel("Repetitions")
plt.plot(aveA,label='A')
plt.plot(aveB,label='B')
plt.plot(ave,label='Both')
plt.legend()
Out[7]:
<matplotlib.legend.Legend at 0x16f3866af60>

Conclusions

As the plot shows that the game a and the game b goes down, it means that are losing games, in contrast to the green one that, on averages, is a wining strategy.

Random Numbers: Accepting or Rejecting

We define the gaussian we want to sample

In [8]:
def gauss(x):
    return np.exp(-0.5*x**2)

We repeat the same procedure until 50000 are acepted

In [9]:
data=np.array([])
while len(data)<50000:
    x=np.random.random()*6-3
    y=np.random.random()
    if y<gauss(x):
        data=np.append(data,x)

we plot the results

In [10]:
hist=plt.hist(data,density=True,bins=100)

We separate the bins and amplitude information

In [11]:
x=(hist[1][:-1]+hist[1][1:])/2
y=hist[0]

We import the fit function and define the model for the fit

In [12]:
from scipy.optimize import curve_fit
def model(x,a):
    return a*gauss(x)

We fit

In [13]:
popt,pcov=curve_fit(model,x,y)

we plot the final results.

In [14]:
plt.hist(data,density=True,bins=100)
plt.plot(x,model(x,*popt))
Out[14]:
[<matplotlib.lines.Line2D at 0x16f38f314e0>]

The percentage error of the fit is

In [15]:
print(pcov.diagonal()**0.5/popt)

[0.00381698]

So, our model seems to be good to generate these kind of points.

Random Numbers on Sphere

If one generates the random numbers (on the apropiate domains),

In [16]:
theta=np.random.random(1000)*np.pi
phi=np.random.random(1000)*2*np.pi
r=1

We get the following plot

In [17]:
plt.scatter(theta,phi)
Out[17]:
<matplotlib.collections.PathCollection at 0x16f390e5d30>

But, it does not represent the sphere, so let us go to the cartesian coordinates

In [18]:
x,y,z=r*np.cos(phi)*np.sin(theta),r*np.sin(phi)*np.sin(theta),r*np.cos(theta)

So when we plot a projection of these points

In [19]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(x,y)
Out[19]:
<matplotlib.collections.PathCollection at 0x16f38f1cef0>

Or other

In [20]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(x,z)
Out[20]:
<matplotlib.collections.PathCollection at 0x16f3a17cc50>

Or other

In [21]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(y,z)
Out[21]:
<matplotlib.collections.PathCollection at 0x16f3a1ddba8>

We do not have the same results, and if they are uniform we shall not expect to have preferred directions,

In [22]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_aspect('equal')
ax.scatter(x,y,z)
Out[22]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x16f3a2853c8>

Gaussian

In [23]:
N=1000
x=np.random.normal(0,1,N)
y=np.random.normal(0,1,N)
z=np.random.normal(0,1,N)
norm=np.sqrt(x**2.+y**2.+z**2.)
x/=norm
y/=norm
z/=norm
In [24]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(x,y)
Out[24]:
<matplotlib.collections.PathCollection at 0x16f3a372898>
In [25]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(x,z)
Out[25]:
<matplotlib.collections.PathCollection at 0x16f3a3d81d0>
In [26]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(y,z)
Out[26]:
<matplotlib.collections.PathCollection at 0x16f3a420860>
In [27]:
fig=plt.figure()
ax=fig.add_subplot(111,projection='3d')
ax.set_aspect("equal")
ax.scatter(x,y,z)
Out[27]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x16f3a485978>

Jacobian

In [28]:
N=1000
phi=np.random.random(N)*2*np.pi
theta=np.arccos(2*np.random.random(N)-1)

r=1

x,y,z=r*np.cos(phi)*np.sin(theta),r*np.sin(phi)*np.sin(theta),r*np.cos(theta)
In [29]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(x,y)
Out[29]:
<matplotlib.collections.PathCollection at 0x16f3a572438>
In [30]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(x,z)
Out[30]:
<matplotlib.collections.PathCollection at 0x16f3a5cdd30>
In [31]:
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.scatter(y,z)
Out[31]:
<matplotlib.collections.PathCollection at 0x16f3a6326a0>
In [32]:
fig=plt.figure()
ax=fig.add_subplot(111,projection='3d')
ax.set_aspect("equal")
ax.scatter(x,y,z)
Out[32]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x16f3a69cc18>

Metropolis Hastings

In [33]:
mu,sig,N = 0,1,1000000

def q(x):
    return (1./(np.sqrt(2*np.pi*sig**2)))*(np.exp(-((x-mu)**2)/(2*sig**2)))

def metropolis(N):
    r = 0
    p = q(r)
    points = []
    
    for i in range(N):
        rn =np.random.uniform(-10,10)
        pn = q(rn)
        
        if pn >= p:
            p = pn
            r = rn
        else:
            u = np.random.rand()
            if u < pn/p:
                p = pn
                r = rn
        points.append(r)    
    points = np.array(points)
    return points
    
    
pts=metropolis(N)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(pts, bins=1000);

Single Random Walk and Rolling Average

In [34]:
x=np.random.normal(0,1,100)

def cumulative_sum(x):
    cum=np.array([x[0]])
    for i in range(1,len(x)):
        cum=np.append(cum,cum[i-1]+x[i])
    return cum
In [35]:
test=cumulative_sum(x)
x=x.cumsum()
In [36]:
print(test-x)

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0.]
In [37]:
plt.plot(x)
Out[37]:
[<matplotlib.lines.Line2D at 0x16f3b0a8fd0>]
In [38]:
lin=np.linspace(0,100,100)
noi=2*np.sin(lin/5)+x
In [39]:
plt.plot(noi)
Out[39]:
[<matplotlib.lines.Line2D at 0x16f3b235a58>]
In [40]:
def roll(x,w):
    rolling=np.array([])
    for i in range(w,len(x)):
        rolling=np.append(rolling,x[i-w:i].mean())
    return rolling
In [41]:
rolling=roll(noi,5)
In [42]:
plt.plot(noi)
plt.plot(np.linspace(5,100,len(rolling)),rolling)
Out[42]:
[<matplotlib.lines.Line2D at 0x16f3b1be128>]

Random Walks: DIffusion

In [43]:
number=10000
steps=1000
x=np.zeros(number)
y=np.zeros(number)
stdx=[]
stdy=[]

for step in range(steps):
    x+=np.random.normal(0,1,number)
    y+=np.random.normal(0,1,number)
    stdx.append(np.var(x)**0.5)
    stdy.append(np.var(y)**0.5)
    
fig=plt.figure()
ax=fig.add_subplot(111)
ax.set_aspect("equal")
#ax.hexbin(x,y,gridsize=30);
ax.hist2d(x,y,bins=np.linspace(-100,100,100));
#plt.savefig("plots/"+str(step)+".png")
#plt.close()
In [44]:
plt.plot(stdx)
plt.plot(stdy)
Out[44]:
[<matplotlib.lines.Line2D at 0x16f3b292470>]

Information Theory: Entropy

In [45]:
from scipy.stats import entropy
print(entropy(hist[0]))

4.219687484693533
In [46]:
def entropy2(P):
    ent=0.
    P=P/np.sum(P)
    for pi in P:
        if pi!=0:
            ent-=pi*np.log(pi)
    return ent
In [ ]:
In [47]:
x=np.random.random(1000000)
hist=plt.hist(x,density=True,bins=100);
print(hist[0].sum())
print(entropy(hist[0]))
print(entropy2(hist[0]))

100.00012015338218
4.605125611116569
4.6051256111165655
In [48]:
x=np.random.normal(0,1,1000000)
hist=plt.hist(x,density=True,bins=100);

print(entropy(hist[0]))
print(entropy2(hist[0]))

3.7332430659669873
3.7332430659669864
In [49]:
x=np.random.poisson(15000,1000000)
hist=plt.hist(x,density=True,bins=100);

print(entropy(hist[0]))
print(entropy2(hist[0]))

3.736089406862346
3.7360894068623454

The distribution with the largest entropy is the uniform distribution