Yes, you can simulate it using a big transition matrix and roll it, while observing the cumulative probability of jumping back to the corner, an absorbing state. This is a Markov traverse on a graph with 64 nodes and a lot of sides, each representing allowed move of a knight on a board.
The average is suprisingly high: 168. This many steps in average is the length of the cycle in this graph when you start in the corner.
Here's the solution in Python:
import numpy as np
def jump(i,j):
"true if jump is allowed"
ret = (abs(i) == 2 and abs(j) == 1) or (abs(i) == 1 and abs(j) == 2)
return ret
# init state after 1st jump
s1 = np.zeros(64) # state after 1st jump
for i in range(8):
for j in range(8):
if jump(i,j):
s1[i*8 + j] = 1
s1 = s1 / np.sum(s1)
# print(np.reshape(s1,(8,8)))
p = np.zeros((64,64)) # transition matrix, where state num k = 8 * i + j
# init tx matrix
for ki in range(8):
for kj in range(8):
k = ki*8 + kj
for i in range(8):
for j in range(8):
if jump(ki-i,kj-j):
p[k, i*8 + j] = 1
p[k, :] = p[k, :] / np.sum(p[k, :])
# print(ki,kj,np.reshape(p[k,:],(8,8)))
p[0,:] = np.zeros(64)
p[0,0] = 1
# print(0,0,np.reshape(p[0,:],(8,8)))
steps = 10000
cumps = np.zeros(steps) # cumulative probability of each cycle length
s = s1
for step in range(2,steps):
s = np.matmul(s, p)
cumps[step] = s[0]
# print(step, ps[step])
ps = np.diff(cumps) # probability of a step
print("check add up to 1: ",cumps[steps-1])
print("mean steps: ",np.dot(range(1,steps),ps))
import matplotlib.pyplot as plt
# plt.plot(range(1,1001), ps[0:1000])
plt.plot(range(1,1001), cumps[1:1001])
plt.xlabel("cycle length")
plt.ylabel("probability")
plt.title("CDF")
plt.show()
Best Answer
Using the Wikipedia notation, the full matrix of transition probabilities is written as $$P = \left(\begin{array}{cc} Q & R \\ 0 & I \end{array}\right)$$ with the states organized so that the last states are all absorbing. In the OP's example from population genetics, the states would be organized as $1, 2, \ldots, M-1, 0, M$.
If $\tau$ denotes the first time the Markov chain enters the set of absorbing states, the absorption probability for state $k$ given that the Markov chain starts in state $i$ is $$B_{ik} := P(X_{\tau} = k \mid X_0 = i).$$ Observe that if $X_n = j$ for a transient state $j$ then $\tau > n$, and by the time homogeneity of the Markov chain $$P(X_{\tau} = k \mid X_n = j) = B_{jk}.$$ From Wikipedia it follows that $$B_{ik} = (NR)_{ik}.$$
The expected number of visits to state $j$ conditionally on the chain starting in $i$ and being absorbed in $k$ is $$\xi_{ik}(j) = E\left( \sum_{n=0}^{\infty} 1(X_n = j) \ \middle| \ X_0 = i, X_{\tau} = k\right) = \sum_{n=0}^{\infty} P(X_n = j \mid X_0 = i, X_{\tau} = k) .$$ Now the probabilities in the infinite sum can be rewritten as follows \begin{align*} P(X_n = j \mid X_0 = i, X_{\tau} = k) & = \frac{P(X_{\tau} = k, X_n = j \mid X_0 = i)}{P(X_{\tau} = k \mid X_0 = i)} \\ & = \frac{P(X_{\tau} = k \mid X_n = j, X_0 = i)P(X_n = j \mid X_0 = i)}{P(X_{\tau} = k \mid X_0 = i)} \\ & = \frac{P(X_{\tau} = k \mid X_n = j)}{P(X_{\tau} = k \mid X_0 = i)} (Q^n)_{ij} \\ & = \frac{B_{jk}}{B_{ik}}(Q^n)_{ij}. \end{align*} From this we obtain the formula \begin{align*} \xi_{ik}(j) & = \frac{B_{jk}}{B_{ik}}\sum_{n=0}^{\infty} (Q^n)_{ij} = \frac{B_{jk}}{B_{ik}}N_{ij}, \end{align*} where the absorption probabilities $B_{ij}$ and $B_{ik}$ can be computed from $N$ and $R$ using the formula above.