それでは、メトロポリスヘイスティングス法を使って、まんべんなく割り当てするパターンを生成するプログラムを作っていきます。
シャフル
まずは、ランダムな初期パターンをつくる関数です。１番目の要素をランダムに選んだ要素と交換し、次に、２番目をランダムな要素と交換、という具合に全ての要素を混ぜ合わせます。
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def suffle(s):
"""
s をシャッフル
"""
s1 = s.copy()
for i in range(len(s)):
j = i
while j==i:
j = np.random.randint(len(s))
wk = s1[i]
s1[i] = s1[j]
s1[j] = wk
return s1
# test -----------------
s_init = np.array([1, 2, 3, 4, 5, 6, 7, 8])
ss = []
for i in range(5):
ss.append(suffle(s_init))
ss = np.array(ss)
print(ss)
output =:
[[5 2 6 7 3 8 1 4]
[1 6 5 7 3 2 8 4]
[2 4 5 6 3 7 8 1]
[5 1 4 2 7 6 8 3]
[8 3 4 5 6 2 1 7]]
テストでは、1から8までのリストの順番をランダムにしています。これを、5回試しています。
状態遷移関数
あるパターンsから次のパターンを生成する関数です。
def next_state(s):
"""
sの二つの要素を入れ替えて
sとは異なる配列にする
"""
s1 = s.copy()
i = np.random.randint(len(s))
j = i
#while s1[j] == s1[i] :
# j = np.random.randint(len(s))
while j == i :
j = np.random.randint(len(s))
wk = s1[i]
s1[i] = s1[j]
s1[j] = wk
return s1
# test ---------------
s_init = np.array([1, 2, 3, 4, 5, 6, 7, 8])
print('init state=',s_init)
print('next state=',next_state(s_init))
output =:
init state= [1 2 3 4 5 6 7 8]
next state= [8 2 3 4 5 6 7 1]
テストでは、1番目と8番目の要素が交換されました。
エネルギー関数
過去のパターンの履歴ssを使って、sのエネルギーを計算する関数です。
def energy(ss, s):
"""
エネルギー関数
"""
c = []
for i in range(ss.shape[1]):
c.append(np.sum(ss[:, i] == s[i]))
c = np.array(c)
eng = np.sum(c)
return eng
# test
ss=np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[1, 0, 0, 0]])
s=np.array([1, 0, 0, 0])
print(ss)
print('---------')
print(s)
print('energy=',energy(ss, s))
output=:
[[1 0 0 0]
[0 1 0 0]
[1 0 0 0]]
———
[1 0 0 0]
energy= 10
エネルギーは、過去のパターンssとの重なりの数ですので、10で間違いないですね。
q関数
履歴ssに対して、パターンsの規格化していない確率qを返す関数です。内部でenergy()が使われています。
def prob_q(s, ss=None, beta=1.0):
if ss is None:
return 1
else:
e = energy(ss, s)
q = 2**(-beta*e) # p1 / p0
return q
# test -----------------
ss = np.array([[1,0,0,0],[0,1,0,0],[1,0,0,0]])
s = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
q = []
for i in range(4):
q.append(prob_q(s[i,:], ss, beta=1.0))
q = np.array(q)
p = q / np.sum(q)
print('q=', np.round(q,4))
print('p=',np.round(p,4))
output =:
q= [0.001 0.0039 0.0156 0.0156]
p= [0.027 0.1081 0.4324 0.4324]
テストでは、1000から、0001までの4つ全てのqとpを計算しました。
ボルツマン分布には、自然対数の定eを使いますが、ここでは、2としています。確かめていませんが、計算が早くなるかなと思ってです。
p1/p0 を求める関数
def p1_over_p0(s0, s1, ss=None, beta=1.0):
if ss is None:
return 1
else:
e0 = energy(ss, s0)
e1 = energy(ss, s1)
p10 = 2**(-beta*(e1-e0)) # p1 / p0
return p10
# test -----------------
ss = np.array([[0 , 1, 1]])
s0 = np.array([1, 0, 0])
s1 = np.array([0 , 0 , 1])
print('energy(s0)=', energy(ss, s0))
print('energy(s1)=', energy(ss, s1))
print('prob(p1)/prob(p0)=', p1_over_p0(s0, s1, ss=ss, beta=1.0))
output=:
energy(s0)= 0
energy(s1)= 2
prob(p1)/prob(p0)= 0.25
エネルギーが低いほど確率は大きくなるので、良さそうですね。
パターン生成関数
以上の関数を使ってメトロポリスヘイスティングス法が実装できます。
def sample_s(ss=None, beta=1.0, N=100):
s0 = suffle(ss[0])
if ss is None:
return s0
else:
cnt = 0
while cnt < N:
s1 = next_state(s0)
p10 = p1_over_p0(s0, s1, ss, beta=beta)
if p10 >= 1:
s_next = s1
elif p10 > np.random.rand():
s_next = s1
else:
s_next = s0
s0 = s_next
cnt += 1
return s0
# test ----------------
ss = np.array([[1, 0, 0, 0]])
s = []
for i in range(10):
s.append(sample_s(ss, beta=2.0, N=10))
s = np.array(s)
print(s)
print('------------')
print('mean=', np.round(np.mean(s==1, axis=0), 2))
output =:
[[0 0 1 0]
[0 0 0 1]
[0 1 0 0]
[0 1 0 0]
[0 0 1 0]
[0 1 0 0]
[0 0 1 0]
[0 0 1 0]
[0 1 0 0]
[0 1 0 0]]
————
mean= [0. 0.5 0.4 0.1]
テストは、
1000の履歴に対してパターンを10回生成しました。履歴は固定しています。1000が出にくくなっていることがわかります。
def sample_ss(init_s, Beta, N):
"""
N回まとめてサンプル
"""
ss = init_s.copy();
for i in range(N-1):
s0 = sample_s(ss, Beta)
s0 = np.array([s0])
ss = np.r_[ss, s0]
return ss
# test ------------------
init_s=np.array([[1, 0, 0, 0]])
N = 8
Beta = 1.0
ss = sample_ss(init_s, Beta, N)
print(ss)
print('------------')
print('beta=', Beta)
print('mean=', np.round(np.mean(ss==1, axis=0), 3))
output =:
[[1 0 0 0]
[0 0 1 0]
[0 0 0 1]
[0 1 0 0]
[0 0 0 1]
[0 0 1 0]
[0 1 0 0]
[0 1 0 0]]
————
beta= 1.0
mean= [0.125 0.375 0.25 0.25 ]
テストは、連続的にパターンを8回生成いた場合です。生成したバターンを履歴にいれながら、次のパターンを生成しています。
完全に振り分けがなされると、１つの列に2つの1が出てきますが、確率的なので、少しずれています。
betaは1.0の場合です。
# test -------------------
init_s=np.array([[1, 0, 0, 0]])
N = 8
Beta = 4.0
ss = sample_ss(init_s, Beta, N)
print(ss)
print('------------')
print('beta=', Beta)
print('mean=', np.round(np.mean(ss==1, axis=0), 3))
output =:
[[1 0 0 0]
[0 0 0 1]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]
[0 1 0 0]
[0 0 1 0]
[1 0 0 0]]
————
beta= 4.0
mean= [0.25 0.25 0.25 0.25]
betaを大きくすると、振り分けをきっちりやるようになるので、全ての列で2回ずつ1が出現しました。
いいかんじですね。
このアルゴリズムをjavascriptで実装してGUIも作ったのが、
でした。