<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Problem 2
Define $h(a) = P(\tau_{\phi} < \tau_{\dagger}) | X_0=a)$
To show
Step1:  T = 0.05
Step2:  T=10
Step3:  Correlation plots
Step4:  Result
The autocorrelation seems to be high even for large values of $N_{step}$ for both the temperature values. I expected higher $T$ to yield lower autocorrelations.
Problem 1
Let the state space be $S = {\phi, \alpha, \beta, \alpha+\beta, pol, \dagger}$ 
Definitions
Step5:  $\alpha$
Step6:  $\beta$
Step7:  $\alpha+\beta$
Step8:  pol
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1)
D = np.random.rand(100,100)
## This is not symmetric, so we make it symmetric
D = (D+D.T)/2
print (D)
import math
N_steps = 10000
def L(sigma):
    s=0
    for i in range(0, len(sigma)-1):
        s+=D[sigma[i], sigma[i+1]]
    return s
def propose(sigma):
    r = np.random.choice(len(sigma), 2)
    rs = np.sort(r)
    j,k=rs[0],rs[1]
    x=(sigma[j:k])#.reverse()
    x=x[::-1]
    x0=  sigma[:j]
    x1 = sigma[k:]
    y=np.concatenate((x0,x,x1))
    return y
def pi(sigma,T):
    return math.exp(-L(sigma)/T)
def metropolis(sigma,T,L_0):
    sigma_n = propose(sigma)
    L_n = L(sigma_n)
    pi_ab = math.exp(-(L_n-L_0)/T)
    q = min(1, pi_ab)
    b = np.random.uniform(size=1)
    if (b<q):
        return sigma_n
    else:
        return sigma
    
sigma_0 = np.random.choice(100,100)
L_0 = L(sigma_0)
print sigma_0
T = [0.05,10]
def plotter(t):
    L_history = []
    sigma_history = []
    sigma_0 = np.random.choice(100,100)
    L_0 = L(sigma_0)
    L_history.append(L_0)
    sigma_history.append(sigma_0)
    sigma = metropolis(sigma_0,t,L_0)
    for i in range(1, N_steps):
        sigma_t = metropolis(sigma_history[i-1],t,L_history[i-1])
        L_1 = L(sigma_t)
        L_history.append(L_1)
        sigma_history.append(sigma_t)
    plt.figure(0)
    plt.hist(L_history, 20)
    #plt.xlim(min(L_history)-25, max(L_history)+0.5)
    plt.xlabel('Length')
    plt.ylabel('Frequency')
    plt.title('Frequency of L')
    plt.figure(1)
    plt.plot(range(1, N_steps+1),L_history)
    plt.ylim(min(L_history), max(L_history))
    plt.xlabel('N_steps')
    plt.ylabel('L')
    plt.title('Variation of L with N_steps')
    return L_history
Explanation: Problem 2
Define $h(a) = P(\tau_{\phi} < \tau_{\dagger}) | X_0=a)$
To show: $P(X_{n+1}=b|X_n=a\ and\ \tau_{\phi} < \tau_{\dagger}) = \frac{h(b)}{h(a)}P_{ab}$
LHS: $P(X_{n+1}=b|X_n=a\ and\ \tau_{\phi} < \tau_{\dagger})$ 
Incorrect version
LHS: 
~~\begin{align}
~~P(X_{1}=b|X_0=a\ and\ \tau_{\phi} < \tau_{\dagger}) &= P(X_{1}=b|X_0=a, \tau_{\phi} < \tau_{\dagger}|X_0=a)~~\
&=P(\tau_{\phi} < \tau_{\dagger}|X_0=a) \times P(X_1=b|X_0=a)\ \text{if $a \notin \phi,\dagger$}\
&=h(a) \times P_{ab} \
&\neq RHS??
\end{align}~~
Corrected Version
Consider 
$P(A|B,C)=P(X_{n+1}=b|X_n=a\ and\ \tau_{\phi} < \tau_{\dagger})$
Then
$P(A|B,C)= \frac{P(A,B,C)}{P(B,C)}=\frac{P(A,C|B)P(B)}{P(C|B)P(B)}=\frac{P(A,C|B)}{P(C|B)} =\frac{P(A|C) \times{P(B|A,C)}}{P(B|C)}$
Thus,
$$
LHS=P(X_{n+1}=b|X_n=a\ and\ \tau_{\phi} < \tau_{\dagger}) = \frac{P(X_{n+1}=b,\tau_{\phi} < \tau_{\dagger}|X_n=a)}
{P(\tau_{\phi} < \tau_{\dagger}|X_n=a)}
$$
Now,  $n > min(\tau_{\phi},\tau_{\dagger})$
and hence:
$$
\begin{align}
P(X_{n+1}=b|X_n=a\ and\ \tau_{\phi} < \tau_{\dagger}) &= \frac{P(X_{n+1}=b,\tau_{\phi} < \tau_{\dagger}|X_n=a)}
{P(\tau_{\phi} < \tau_{\dagger}|X_n=a)} \
&=  \frac{P(X_{n+1}=b|X_n=a)\times P(\tau_{\phi} < \tau_{\dagger}|X_n=a,X_{n+1}=b)}{P(\tau_{\phi} < \tau_{\dagger}|X_n=a)}
\end{align}
$$
Using markov property and time homogeneity:
$P(\tau_{\phi} < \tau_{\dagger}|X_n=a,X_{n+1}=b)=P(\tau_{\phi} < \tau_{\dagger}|X_0=b)$
and hence:
$$
\begin{align}
P(X_{n+1}=b|X_n=a\ and\ \tau_{\phi} < \tau_{\dagger}) &=  \frac{P(X_{n+1}=b|X_n=a)\times P(\tau_{\phi} < \tau_{\dagger}|X_n=a,X_{n+1}=b)}{P(\tau_{\phi} < \tau_{\dagger}|X_n=a)}\
&=\frac{h(b)\times P_{ab}}{h(a)}\
&= RHS
\end{align}
$$
Problem 3
If initial state $X_t=A$,
$P(X_{t+1}=A) = 0.5$ and $P(X_{t+1}=A\cup{b}-{a})=0.5$
Observation 1: $X_t$ is irreducbile. The construction allows to reach every state from any state.
Example: Let $A={1,2,3,4,5}$ for $n=10$ and $k=5$. let $a$=3 and let $b=6$
Then we have: $P(X_{t+1}={1,2,3,4,5}) = 0.5$ and $P(X_{t+1}={1,2,4,5,x})=0.5*1/5$ where $x\ \in {6,7,8,9,10}$
Observation 2: For $X$ to be aperioidic, it is imporatant to have the $X_{t+1}=X_t$ with probabulity 0.5(any non-zero probability would do). Otherwise the diagonal of the trasition probability matrix will be zero, and in such cases it is possible for the chain to be periodic. An example (without taking into account the actual transition probabilities) is:
For state space.${1,2,3,4}$
$$
P = \begin{bmatrix}
0 & 0.5 & 0 & 0.5\
0.5 & 0 & 0.5 & 0\
0 & 0.5 & 0 & 0.5\
0.5 & 0 & 0.5 & 0\
\end{bmatrix}
$$
and $(P^2)_{ii}>0$
It is possible to return to the same state with a period of 2:
$P(X_n=2|X_0=1)= 0 \ for\ $n=2k$\ and\ 1\ for\ $n=2k-1$\ where\ k=1,2,3...$
About uniform stationary distribution
From observations 1,2 we know that the markv chain is irreducible and aperiodic. There is another observation:
Observation 3: $P$ the transition probabilty matrix is symmetric.
$P_{ii} = 0.5$
$P_{ij} = 0.5 * \underbrace{\frac{1}{|k|}}\text{Probability of selecting 'i' uniformly}  * \underbrace{\frac{1}{|A|-|k|}}\text{Probability of selecting 'j' uniformly}$ $\forall j \neq i$
and hence $P_{ij} =P_{ji}$ $\implies$ $P=P^T$ $\implies$ $\pi$ is uniformly distributed (Because $P$ is reversible)
Problem 4
Part (4a)
End of explanation
L_t0=plotter(T[0])
Explanation: T = 0.05
End of explanation
L_t1= plotter(T[1])
Explanation: T=10
End of explanation
from scipy.signal import correlate
def autocorr(x):
    xunbiased = x-np.mean(x)
    xnorm = np.sum(xunbiased**2)
    acor = np.correlate(xunbiased, xunbiased, "same")/xnorm
    #result = correlate(x, x, mode='full')
    #result /= result[result.argmax()]
    acor = acor[len(acor)/2:]
    return acor#result[result.size/2:]
cov_t0 = autocorr(L_t0)
cov_t1 = autocorr(L_t1)
plt.plot(cov_t0)
plt.ylabel('Autocorrelation')
plt.xlabel('N_steps')
plt.title('Autocorrelation of L_i for T=0.05')
plt.plot(cov_t1)
plt.ylabel('Autocorrelation')
plt.xlabel('N_steps')
plt.title('Autocorrelation of L_i for T=10')
Explanation: Correlation plots
End of explanation
k_a=0.2
k_b=0.2
k_p=0.5
P = np.matrix([[1-k_a-k_b, k_a ,k_b, 0, 0, 0],
               [k_a, 1-k_a-k_b, 0, k_b, 0, 0],
               [k_b, 0, 1-k_a-k_b, k_a, 0, 0],
               [0, k_b, k_a, 1-k_a-k_b-k_p, k_p, 0],
               [0, 0, 0, 0, 0, 1],
               [0, 0, 0, 1, 0, 0]])
Q=P[1:5,1:5]
iq = np.eye(4)-Q
iqi = np.linalg.inv(iq)
print(iq)
print(iqi)
print 'U={}'.format(iqi[:,-1])
u=iqi[:,-1]
PP = {}
states = ['phi', 'alpha', 'beta', 'ab', 'pol', 'd']
PP['phi']= [1-k_a-k_b, k_a ,k_b, 0, 0, 0]
PP['alpha'] = [k_a, 1-k_a-k_b, 0, k_b, 0, 0]
PP['beta'] = [k_b, 0, 1-k_a-k_b, k_a, 0, 0]
PP['ab']= [0, k_b, k_a, 1-k_a-k_b-k_p, k_p, 0]
PP['pol']=  [0, 0, 0, 0, 0, 1]
PP['d']= [0, 0, 0, 1, 0, 0]
def h(x):
    s=0
    ht=0
    cc=0
    for j in range(1,100):
        new_state=x
        for i in range(1,10000):
            old_state=new_state
            probs = PP[old_state]
            z=np.random.choice(6, 1, p=probs)
            new_state = states[z[0]]
            s+=z[0]
            if new_state=='d':
                ht+=i
                cc+=1
                break
            else:
                continue
    return s/1000, ht/cc
Explanation: Result
The autocorrelation seems to be high even for large values of $N_{step}$ for both the temperature values. I expected higher $T$ to yield lower autocorrelations.
Problem 1
Let the state space be $S = {\phi, \alpha, \beta, \alpha+\beta, pol, \dagger}$ 
Definitions:
1. $\tau_a = { n \geq 0: X_n=a}$
$N = \sum_{k=0}^{\tau_{\phi}}I_{X_k=\dagger}$ 
$u(a) = E[N|X_0=a] \forall a \in S $
$u(a) = \sum_{k=0}^{\tau_{\phi}}P(X_k=\dagger|X_0=a)=\sum_{b \neq a, \dagger }P(X_1=b|X_0=a)P(X_k=\dagger|X_0=b)$ $\implies$ $u(a)=\sum_{b \neq a, \dagger} P_{ab}u(b)$
And hence $u$ solves the following set of equations:
$u=(I-P_{-})^{-1}v$ where v is (0,0,0,1) in this case. and $P_{-}$ represents the matrix with that last and first row and columns removed.
End of explanation
print('Simulation: {}\t Calculation: {}'.format(h('alpha')[1],u[0]))
Explanation: $\alpha$
End of explanation
print('Simulation: {}\t Calculation: {}'.format(h('beta')[1],u[1]))
Explanation: $\beta$
End of explanation
print('Simulation: {}\t Calculation: {}'.format(h('ab')[1],u[2]))
Explanation: $\alpha+\beta$
End of explanation
print('Simulation: {}\t Calculation: {}'.format(h('pol')[1],u[3]))
Explanation: pol
End of explanation
<END_TASK>