Markov Chains and Monte Carlo Methods Computer Practical 1. Markov chains

Ioana A. Cosma and Ludger Evers, Markov Chains and Monte Carlo Methods
http://users.aims.ac.za/~ioana/notes.pdf
CC-by-nc-sa 2.5
http://creativecommons.org/licenses/by-nc-sa/2.5/za/legalcode
Q and A

task1

def m_trans(K, m):
    # A, Bがnp.array型のときA * BやA**nはアダマール積として扱われるので,
    # np.matrix型に変換する.
    K = np.matrix(K)
    return K**m

m_trans(np.array([[0.9, 0.1], [0.3, 0.7]]), 3)

->

matrix([[ 0.804,  0.196],
        [ 0.588,  0.412]])

task2

def dist_t(K, t, lambda_0):
    return np.transpose(lambda_0).dot(m_trans(K, t))

# phone lineの例で,stationary distributionを試すと
dist_t(np.array([[0.9, 0.1], [0.3, 0.7]]), 100, np.array([0.75, 0.25]))

->

matrix([[ 0.75,  0.25]])

task3

def stationary(K):
    K = np.matrix(K)
    dim = len(K)
    b = np.zeros(dim+1)
    b[dim] = 1

    up = np.transpose(K - np.matrix(np.identity(dim)))
    mu = np.linalg.lstsq(np.array(np.concatenate((up, [np.ones(dim)]))), b)
    print(mu[0])


stationary(K)

->

[ 0.75  0.25]

task4

t = np.linspace(0, 25)
a = np.random.uniform()
init = np.array([a, 1-a])
def pow_generator(K): # Kのべき乗を生成するジェネレータ
    K = np.matrix(K)
    tmp = np.identity(len(K))
    while True:
        tmp = K.dot(tmp)
        yield tmp

gen = pow_generator(K)
process = [np.array(init.dot(next(gen)))[0] for t in range(25)]

plt.plot(process, marker=".")
plt.ylim([0, 1])
plt.show()

->

task5

def sample_path(K, init, tau):
    states = list(range(len(init))) # ここではstateは0から番号を付ける
    x0 = np.random.choice(states, p=init)
    path = [x0]
    tmp = x0

    for i in range(tau):
        suc = np.random.choice(states, p=K[tmp])
        tmp = suc
        path.append(suc)

    return path

plt.plot(sample_path(K, [0.5, 0.5], 100))
plt.show()

->

task6

path = sample_path(K, [0.5,0.5],20000)
avg = []
for i in range(1, 100):
    avg.append(np.average(path[: 200 * i]))
plt.plot(avg)
plt.show()

->