#!/usr/bin/env python
# coding: utf-8

# # [Homework Week 8](https://tsoo-math.github.io/ucl2/2021-HW-week8.html)
# 
# For complete written solutions, please [see](https://tsoo-math.github.io/ucl/QHW3-sols.html).  We will mainly be concerned with the Python code here.

# ## A Markov chain
# 
# 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$.

# In[1]:


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  


# In[2]:


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])


# ## Gibbs sampler
# 
# 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.

# In[3]:


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)


# In[ ]:





# ## Simple card shuffling

# In[6]:


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))
    
        
    
    


# ## Endnotes
# 
# Use the ipynb [source](https://tsoo-math.github.io/ucl2/2021-hw-week8.ipynb) for the most update version.

# In[5]:


from datetime import datetime
print(datetime.now())


# In[ ]:





# In[ ]:





# In[ ]: