For complete written solutions, please see. We will mainly be concerned with the Python code here.
The following codes a Markov chain in Python; we run the chain for $n$ steps, and when $n$ is large, we find that the total number of visits $v(i)$ to a state $i$ is such that $v(i)/n \approx \pi(i)$, where $\pi$ is the stationary distribution, which agrees with a version of law of large number for Markov chains. We also check a version of the law of large numbers for the number of transitions from state $1$ to $3$.
import numpy as np
P = np.array( [np.array([1/4, 1/4, 1/2, 0,0]),
np.array([1/4, 1/4, 0,0,1/2]), np.array([1/4,0,0,1/2,1/4]), np.array([0,0,0,1/2, 1/2]), np.array([1/3,1/3,0,0,1/3] ) ] )
print(P)
print(P[1])
print(P[4])
def step(i): # advancing the MC by one step, when you are at state i
q = P[i-1] # this the transition vector for state i
x=-1
u = np.random.uniform() # imagine the interval [0,1] split up into smaller intervals with probabilities q_j summing to 1
j=0
cumq = np.cumsum(q)
while(x==-1):
j = j+1
if (u <= cumq[j-1]): # if u lands in the interval of length q_j, then we jump to state j
x = j
return x
def steps(i,n): # advances the MC by n steps, starting at state i, keeping a complete history
x = np.array([i])
for j in range(n):
x = np.append(x, step(x[j]) )
return x
# print(steps(1,100))
mc = steps(1,100000)
freq = np.unique(mc,return_counts = True) # returns a frequency count
print(freq)
stat = (1/100000)*np.array( [ freq[1][0] , freq[1][1] , freq[1][2] , freq[1][3] , freq[1][4] ]) # approximation of the stationary distribution
print(stat)
stat
np.matmul(stat,P) - stat # checking
[[0.25 0.25 0.5 0. 0. ] [0.25 0.25 0. 0. 0.5 ] [0.25 0. 0. 0.5 0.25 ] [0. 0. 0. 0.5 0.5 ] [0.33333333 0.33333333 0. 0. 0.33333333]] [0.25 0.25 0. 0. 0.5 ] [0.33333333 0.33333333 0. 0. 0.33333333] (array([1, 2, 3, 4, 5]), array([24435, 21391, 12233, 12244, 29698], dtype=int64)) [0.24435 0.21391 0.12233 0.12244 0.29698]
array([-2.09166667e-04, -3.51666667e-04, -1.55000000e-04, -5.50000000e-05, 7.70833333e-04])
def onethree(n):
mc= steps(1,n)
onethree=0
for i in range(n):
if (mc[i]==1 and mc[i+1]==3):
onethree = onethree +1
return onethree
print(onethree(50000)/50000 - P[1-1, 3-1]*stat[1-1])
0.000784999999999994
Our analysis here will be slightly different than in the R version; instead of running the Gibbs sampler for a $100$ steps, recording the output, which should be close to being from the stationary distribution, and then repeating this procedure $10000$ times, we will just run the Gibbs sampler for $10000$ steps, and appeal to a version of the law of large numbers of Markov chains to recover the stationary distribution. Notice the subtle difference, in the first approach, we appeal to the convergence to the stationary distribution, and by doing repeated independent experiments, appeal to the usual law of large numbers; in the second approach we appeal to the law of large numbers for Markov chains.
def fwz(z):
if z==0:
w = np.random.binomial(1, 2/5,1)
else:
w = np.random.binomial(1, 2/3,1)
return w
def fzw(w):
if w==0:
z = np.random.binomial(1,1/4,1)
else:
z = np.random.binomial(1,1/2,1)
return z
def gibbs(n):
w=np.array([1])
x=w
y = fzw(x[0])
for i in range(n):
x = np.append(x, fwz(y[i]) )
y = np.append(y, fzw(x[i+1]) )
return [x,y]
g=gibbs(10000)
counts11 = 0
for i in range(10000):
if (g[0][i]==1 and g[1][i]==1 ):
counts11 = counts11 +1
print(counts11/10000 -1/4)
counts10 = 0
for i in range(10000):
if (g[0][i]==1 and g[1][i]==0 ):
counts10 = counts10 +1
print(counts10/10000-1/4)
counts00 = 0
for i in range(10000):
if (g[0][i]==0 and g[1][i]==0 ):
counts00 = counts00 +1
print(counts00/10000-3/8)
counts01 = 0
for i in range(10000):
if (g[0][i]==0 and g[1][i]==1 ):
counts01 = counts01 +1
print(counts01/10000-1/8)
0.01090000000000002 -0.0030999999999999917 -0.00019999999999997797 -0.007599999999999996
import random
deck = list( range(1,5))
print(deck) # the deck in order
print(random.sample(range(1,5), 4)) # a random deck of 4 cards
def shuffle(x):
if np.random.binomial(1,1/2,1)==1:
t=random.sample(range(1,5), 2)
a= t[0]
b= t[1]
da = x[a-1]
db = x[b-1]
x[a-1]=db
x[b-1]=da
return x
def coupled(x):
deckx = np.array(x)
decky = np.array(random.sample(range(1,5), 4))
n=0
while( np.array_equal(deckx, decky)==False):
deckx = shuffle(deckx)
dekcy = shuffle(decky)
n = n+1
return n
repeat = [coupled(deck) for _ in range(10000) ]
print(np.mean(repeat))
[1, 2, 3, 4] [2, 1, 3, 4] 32.7245
from datetime import datetime
print(datetime.now())
2021-10-15 00:57:51.429038