## Problem 1




The transition matrix is given by:
$$
\begin{bmatrix}
1-\alpha & \alpha\\
\beta & 1-\beta
\end{bmatrix}
$$

### Part (a)

$\eta = min(n>0, X_n=1)$ given $X_0=0$

*To Prove* $\eta \sim Geom(\alpha)$


$\eta = P(X_0=0,X_1=0, \dots X_{n-1}=1, X_n=1)$

Using the Markov property this can be written as:

$$
\eta = P(X_0=0)P(X_1=0|X_0=0)P(X_2=0|X_1=0)P(X_3=0|X_2=0) \dots P(X_{n-1}=0|X_{n-2}=0)P(X_{n}=1|X_{n-1}=0)
$$

And being time-homogenous, this simplifies to:

$$
\eta = P(X_0=0)\big(P(X_1=0|X_0)\big)^{n-1}\times P(X_1=1|X_0=0)
$$

$\implies$ 

$$
\eta = P(X_0=0)\big(1-\alpha)^{n-1}\alpha = \big(1-\alpha)^{n-1}\alpha
$$

And hence $\eta \sim Geom(\alpha)$


### Part (b)

Spectral decomposition of $P$ and value for $P(X_n=1|X_0=0)$

Spectral decomposition of $P$:

$$
det\begin{bmatrix}
\alpha-\lambda & 1-\alpha\\
1-\beta & \beta-\lambda
\end{bmatrix} = 0
$$

$$
\lambda^2 +(\alpha + \beta-2) \lambda + (1-\alpha -\beta) = 0
$$

Thus, $\lambda_1 = 1$ and $\lambda_2 = 1-\alpha-\beta$

Eigenvectors are given by:

$v_1^T = \big( x_1\ x_1 \big)\ \forall\ x_1 \in R$

and for $\lambda_2$ , $v_2 = \big( x_1\ \frac{-\beta x_1}{\alpha} \big)$

Now using Markov property: $P(X_n=1|X_0=0) = (P^n)_{01}$

Now, 

$P^n = VD^nV^{-1}$

where:
$$
V = \begin{bmatrix}
1 & 1\\
1 & \frac{-\beta}{\alpha}
\end{bmatrix}
$$

and 

$$
D = \begin{bmatrix}
1 & 0 \\
0 & (1-\alpha-\beta)
\end{bmatrix}
$$

$$
V^{-1} = \frac{-1}{\frac{\beta}{\alpha}+1}\begin{bmatrix}
-\frac{\beta}{\alpha} & -1 \\
-1 & 1
\end{bmatrix}
$$

Thus,

$$
P^n = \begin{bmatrix}
1 & 1\\
1 & \frac{-\beta}{\alpha}
\end{bmatrix} \times \begin{bmatrix}
1 & 0 \\
0 & (1-\alpha-\beta)^n
\end{bmatrix} \times \frac{-1}{\frac{\beta}{\alpha}+1}\begin{bmatrix}
-\frac{\beta}{\alpha} & -1 \\
-1 & 1
\end{bmatrix}
$$

$$
P^n = \frac{1}{\alpha+\beta} \begin{bmatrix}
\beta + \alpha(1-\alpha-\beta)^n & \alpha-\alpha(1-\alpha-\beta)^n\\
\beta - \beta(1-\alpha-\beta)^n & \alpha + \beta(1-\alpha-\beta)^n
\end{bmatrix}
$$

## Part (c)

When $\alpha+\beta=1$, the eigen values are $\lambda_1=1$ and $\lambda_2=0$ and hence

$$
P^n = \begin{bmatrix}
\beta & \alpha \\
\beta & \alpha
\end{bmatrix}
$$

*Check:*


Also consider the following identifiy: $P^{n+1}=PP^n$

then:

$$
\begin{bmatrix}
p_{00}^{n+1} & p_{01}^{n+1}\\
p_{10}^{n+1} & p_{11}^{n+1}\\
\end{bmatrix} = \begin{bmatrix} 
p_{00}^n & p_{01}^n\\
p_{10}^n & p_{11}^n
\end{bmatrix} \times \begin{bmatrix}
1-\alpha & \alpha\\
\beta & 1-\beta
\end{bmatrix}
$$

$\implies$

$$
\begin{align}
p_{11}^{n+1} &= p_{10}^n(\alpha) + p_{11}^n(1-\beta)\\
 &= (1-p_{11}^n)(\alpha) +(p_{11}^n)(1-\beta)\\
 &= \alpha + (1-\alpha-\beta)p_{11}^n
\end{align}
$$

Consider the recurrence:

$$
x_{n+1} = \alpha+(1-\alpha-\beta)x_n
$$

Constant solution $x_n=x_{n+1}=x$ is given by: $x=\frac{\alpha}{\alpha+\beta}$

Now let $y_n = x_n-x=x_n-\frac{\alpha}{\alpha+\beta}$ then,

$y_{n+1} = (1-\alpha-\beta)y_n$ and hence $y_n=(1-\alpha-\beta)^n y_0$

Thus,
$$p_{11}^{n} = (1-\alpha-\beta)^np_{11}^0 +\frac{\alpha}{\alpha+\beta}$$

Given $P_{00}=\frac{\beta}{\alpha+\beta}$ and $\alpha+\beta=1$ and hence:

$p_{11}^n = \frac{\alpha}{\alpha+\beta} = \alpha$

and hence, $p_{10}^n = \beta$

Similary,

$p_{00}^n = \beta$ and $p_{01}^n = \alpha$


## Problem 2

$P(X_1=0) \frac{\beta}{\alpha+\beta}$ and hence $P(X_1=1) = \frac{\alpha}{\alpha+\beta}$

$X=X_1X_2\dots X_n$ and $Y=Y_1Y_2\dots Y_n$ representes the reverse string $Y_k=X_{n+k-1}$


### Part (a)

Given string of digits: $a_1,a_2,a_3 \dots a_n $ to find: $P(Y_1=a_1,Y_2=a_2,Y_3=a_3\dots Y_n=a_n)$

$$
\begin{align}
P(Y_1=a_1,Y_2=a_2,Y_3=a_3\dots Y_n=a_n) &= P(X_1=a_n,X_2=a_{n-1}, \dots X_n=a_1) \\
&= P(X_1=a_n)P(X_2=a_{n-1}|X_1=a_n)P(X_3=a_{n-2}|X_2=a_{n-1})\dots P(X_n=a_1|X_{n-1}=a_2) \\
&= P(X_1=a_n)(P_{a_n a_{n-1}})(P_{a_{n-1} a_{n-2}}) \dots (P_{a_2 a_1})
\end{align}
$$

The problem asked about not using spectral decomposition, but I was not sure how spectral decomposition would have come in handy if the states $a_i$ are not specified explicitly.



### Part (b)

$$
z=\begin{cases}
X & if \theta = H\\
Y & otherwise
\end{cases}
$$

Given function f such that, $f :\{0,1\}^n \longrightarrow \{H,T\}$
To show: $P(f(Z)=\theta)=0.5$

$P(\theta=H) = P(\theta=T) = 0.5$

Given Z, guess $\theta$: 

$P(\theta=H|Z=X) = \frac{P(\theta=H, Z=X)}{P(Z=X)}$

Z, has only two possible values: $H$ and $T$ and hence assuming the guess function is unbiased:

$P(f(Z)=H) = P(f(Z)=T)=0.5$


## Problem 3

$$
\tau = min\{ n \geq 0: X_n=\dagger\}
$$

$$
E[\tau] = E_a[E_a[\tau|X_n=a]\ where\ a \in \{\phi, \alpha, \beta, \alpha+\beta, pol, \dagger\}
$$

Let $S=\{\phi, \alpha, \beta, \alpha+\beta, pol, \dagger\}$


Consider for $a\neq \dagger$:
$$
h(a) = E[\tau|X_0=a] = \sum_{s \in S}P_{as} \times (1) + P_{as}\times E[\tau|X_0=s) 
$$

$\implies$ 

$$
h(a) = ((I-P_{-})^{-1})_a
$$

where $P_{-}$ represents the matrix with the row and column representng $X_i=\dagger$ removed.





```python
%matplotlib inline
from __future__ import division
import numpy as np
from numpy import linalg as LA
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 = [[k_a-k_b,k_a,k_b,0,0],
     [k_a,k_a+k_b,0,k_b,0],
     [k_b,0,k_a+k_b,0,0],
     [0,k_b,k_a,k_a+k_b+k_p,k_p],
     [0,0,0,0,0]]
qq = np.array(q)

print(P)
states = ['phi', 'alpha', 'beta', 'ab', 'pol', 'd']

```

    [[ 0.6  0.2  0.2  0.   0.   0. ]
     [ 0.2  0.6  0.   0.2  0.   0. ]
     [ 0.2  0.   0.6  0.2  0.   0. ]
     [ 0.   0.2  0.2  0.1  0.5  0. ]
     [ 0.   0.   0.   0.   0.   1. ]
     [ 0.   0.   0.   1.   0.   0. ]]



```python
import networkx as nx

G=nx.from_numpy_matrix(P,create_using=nx.MultiDiGraph())
G.edges(data=True)
#nx.draw_graphviz(G)# labels=states)
nx.write_dot(G,'G.dot')
```


```python
!neato -T png G.dot > multi.png
```



The markov chain seems to be irreducible
One way to obtain the stationary state is to look at the eigen vectors correspendoing to the eigen value of 1. However, the eigen vectors come out to be imaginary. This seemed to be an issue wwith the solver so I relied on solving the system of equation: $\pi = P\pi$


```python
w, v = LA.eig(P)
for i in range(0,6):
    print 'Eigen value: {}\n Eigen vector: {}\n'.format(w[i],v[:,i])
```

    Eigen value: (-0.382094716832+0.673687371073j)
     Eigen vector: [[ 0.01610056-0.01846064j]
     [-0.00843895+0.07244209j]
     [-0.00843895+0.07244209j]
     [-0.21867795-0.36569039j]
     [ 0.71032012+0.j        ]
     [-0.27140957+0.4785337j ]]
    
    Eigen value: (-0.382094716832-0.673687371073j)
     Eigen vector: [[ 0.01610056+0.01846064j]
     [-0.00843895-0.07244209j]
     [-0.00843895-0.07244209j]
     [-0.21867795+0.36569039j]
     [ 0.71032012-0.j        ]
     [-0.27140957-0.4785337j ]]
    
    Eigen value: (1+0j)
     Eigen vector: [[ 0.40824829+0.j]
     [ 0.40824829+0.j]
     [ 0.40824829+0.j]
     [ 0.40824829+0.j]
     [ 0.40824829+0.j]
     [ 0.40824829+0.j]]
    
    Eigen value: (0.755104836853+0j)
     Eigen vector: [[-0.49079708+0.j]
     [-0.19031250+0.j]
     [-0.19031250+0.j]
     [ 0.34320513+0.j]
     [ 0.60192068+0.j]
     [ 0.45451322+0.j]]
    
    Eigen value: (0.309084596812+0j)
     Eigen vector: [[-0.63710764+0.j]
     [ 0.46336106+0.j]
     [ 0.46336106+0.j]
     [-0.03688671+0.j]
     [-0.38611374+0.j]
     [-0.11934181+0.j]]
    
    Eigen value: (0.6+0j)
     Eigen vector: [[  3.12250226e-17+0.j]
     [  7.07106781e-01+0.j]
     [ -7.07106781e-01+0.j]
     [ -3.12250226e-17+0.j]
     [ -8.67361738e-17+0.j]
     [ -5.20417043e-17+0.j]]
    



```python
## Solve for (I-Q)^{-1}
iq = np.linalg.inv(np.eye(5)-qq)
iq_phi = iq[0,0]
iq_alpha = iq[1,1]
iq_beta = iq[2,2]
iq_alphabeta = iq[3,3]
iq_pol = iq[4,4]

```

EDIT: I made correction to solve for corrected $\pi$, by acounting for $P^T$ and not $P$


```python
A = np.eye(6)-P.T
A[-1,:] = [1,1,1,1,1,1]

B = [0,0,0,0,0,1]
X=np.linalg.solve(A,B)
print(X)

```

    [ 0.2  0.2  0.2  0.2  0.1  0.1]


Stationary state is given by $\pi = (0.1667, 0.1667, 0.1667, 0.1667, 0.1667, 0.1667)$ The mean number of visits per unit time to $\dagger$ are $\frac{1}{\pi_6} = 6$ However strangely this does not satisfy $\pi=P\pi$. I was not able to figure out where I went wrong.

EDIT: I made correction to solve for corrected $\pi$, by acounting for $P^T$ and not $P$, so this no longer holds


```python
#EDIT: I made correction to solve for corrected $\pi$, by acounting for $P^T$ and not $P$
print('\pi*P={}\n'.format(X*P))
print('But \pi={}'.format(X)) 
```

    \pi*P=[[ 0.2  0.2  0.2  0.2  0.1  0.1]]
    
    But \pi=[ 0.2  0.2  0.2  0.2  0.1  0.1]




Simulating the chain:

General strategy: Generate a random number $\longrightarrow$ Select a state $\longrightarrow$ Jump to state $\longrightarrow$ Repeat



```python
## phi
np.random.seed(1)
```


```python
PP = {}
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]
```


```python
##For $h(\phi)$
x0='phi'
x='phi'
def h(x):
    s=0
    new_state=x
    for i in range(1,1000):
        old_state=new_state
        probs = PP[old_state]
        z=np.random.choice(6, 1, p=probs)
        new_state = states[z[0]]
        #print('{} --> {}'.format(old_state, new_state))
        s+=z[0]
    return s/1000

```

## Part (a,b,c)



```python
print(r'$h(\phi)$: From simulation: {}; From calculation: {}'.format(h('phi'),iq_phi))

```

    $h(\phi)$: From simulation: 2.199; From calculation: 1.66666666667



```python
print(r'$h(\alpha)$: From simulation: {}; From calculation: {}'.format(h('alpha'),iq_alpha))

```

    $h(\alpha)$: From simulation: 2.364; From calculation: 7.77777777778



```python
print(r'$h(\beta)$: From simulation: {}; From calculation: {}'.format(h('beta'),iq_beta))

```

    $h(\beta)$: From simulation: 1.988; From calculation: 2.22222222222



```python
print(r'$h(\alpha+\beta)$: From simulation: {}; From calculation: {}'.format(h('ab'),iq_alphabeta))

```

    $h(\alpha+\beta)$: From simulation: 2.189; From calculation: 43.3333333333



```python
print(r'$h(\pol)$: From simulation: {}; From calculation: {}'.format(h('pol'),iq_pol))

```

    $h(\pol)$: From simulation: 2.031; From calculation: 1.0



```python
old_state = [0.1,0.2,0.3,0.4,0,0]

def perturb(old_state):
    new_state = old_state*P
    return new_state
new_state = [0,0,0,0,0,1]

while not np.allclose(old_state, new_state):
    old_state, new_state = new_state, perturb(old_state)
    
print old_state

```

    [[ 0.19999998  0.20000004  0.20000004  0.20000002  0.1000002   0.09999973]]



```python
# EDIT: I made correction to solve for corrected $\pi$, by acounting for $P^T$ and not $P$
print('From calculation(which is NO LONGER wrong!), stationary distribution:{}'.format(X))

print('From simulation, stationary distribution: {}'.format(old_state))
```

    From calculation(which is NO LONGER wrong!), stationary distribution:[ 0.2  0.2  0.2  0.2  0.1  0.1]
    From simulation, stationary distribution: [[ 0.19999998  0.20000004  0.20000004  0.20000002  0.1000002   0.09999973]]



```python

```

0259011da13d6e0d116b37fdf3c120ba43ead78d
ipynb
Jupyter Notebook
2015_Fall/MATH-578B/Homework1/Homework1.ipynb
NeveIsa/hatex
c5cfa2410d47c7e43a476a8c8a9795182fe8f836
MIT
2015-09-10T02:45:33.000Z
2022-02-10T03:20:47.000Z
2015_Fall/MATH-578B/Homework1/Homework1.ipynb
NeveIsa/hatex
c5cfa2410d47c7e43a476a8c8a9795182fe8f836
MIT
2015-09-16T23:11:00.000Z
2015-09-23T21:21:52.000Z
2015_Fall/MATH-578B/Homework1/Homework1.ipynb
saketkc/hatex
c5cfa2410d47c7e43a476a8c8a9795182fe8f836
MIT
2015-09-25T19:06:45.000Z
2022-02-10T03:21:09.000Z
Qwen/Qwen-72B
1. YES
2. YES