怪しいダイス(hmmlearn)

ダイスの出目を隠れマルコフモデルでふにふにすることについて質問されたので書いた。 とりあえずhmmlearnで簡単に。

gistのほうが見やすいかも DirtyDice.ipynb

!pip install hmmlearn
Collecting hmmlearn
  Downloading hmmlearn-0.2.0.tar.gz (107kB)
[K    100% |████████████████████████████████| 112kB 1.8MB/s ta 0:00:01
[?25hBuilding wheels for collected packages: hmmlearn
  Running setup.py bdist_wheel for hmmlearn ... [?25ldone
[?25h  Stored in directory: /home/jovyan/.cache/pip/wheels/de/6c/25/c5fc03ff58b509d7b9769a0d5e2f69fc3fe89feebd70b89b86
Successfully built hmmlearn
Installing collected packages: hmmlearn
Successfully installed hmmlearn-0.2.0
import math
import numpy as np
from hmmlearn import hmm


np.random.seed(1) # デバッグのためにシードを設定しておく
# まずマトモなダイスの出力を再現してみる
n_dice = 1 #ダイスの数は1で
n_dice_faces = 6 # 6面とする
model = hmm.MultinomialHMM(n_components=n_dice)
model.startprob_ = np.array([1.0]) # ダイスは1つなので、最初のダイスが選ばれる確率は1
model.transmat_ = np.array([[1.0]]) # ダイスは1つなので、最初のダイスから常に最初のダイスが選ばれる
model.emissionprob_ = np.array([[1/6, 1/6, 1/6, 1/6, 1/6, 1/6]]) # ダイスのそれぞれの面は等しい確率で選ばれる
# ダイスを100回振ったとして
n_try = 100
x, z = model.sample(n_try) 
x = x.reshape(1, -1)[0] # 見やすく変形しているだけ
# ダイスを振った記録
x
array([5, 4, 1, 2, 2, 1, 2, 4, 2, 4, 5, 1, 5, 4, 4, 4, 2, 0, 0, 1, 5, 3, 2,
       1, 1, 5, 0, 2, 1, 5, 3, 0, 2, 3, 3, 4, 3, 0, 5, 3, 5, 5, 3, 3, 0, 4,
       1, 2, 2, 5, 0, 3, 3, 3, 0, 0, 4, 2, 4, 1, 3, 1, 1, 2, 3, 4, 0, 5, 4,
       3, 4, 4, 5, 1, 2, 3, 4, 2, 2, 3, 0, 0, 0, 4, 4, 2, 2, 2, 0, 0, 0, 0,
       5, 5, 0, 3, 5, 4, 2, 1])
# 常に最初のダイスが選ばれる
z
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0])
# この試行結果からそれぞれの面の確率を求めると、
remodel = hmm.MultinomialHMM(n_components=n_dice, n_iter=100)
remodel.fit(np.array([x]).T)
remodel.emissionprob_
array([[ 0.165,  0.17 ,  0.17 ,  0.165,  0.165,  0.165]])
# さて、ここで、"怪しい" サイコロのログを調べてみるとする  
# Ref: http://www.fward.net/archives/1548
x = "12442552156235244611513665441521334522422516131461324351224626256625324354313241416144163515363423135135313531535155116446122436632362242642626642654256461421553232365354361461451264546335353161532646"
x = [int(s)-1 for s in x]

remodel = hmm.MultinomialHMM(n_components=n_dice, n_iter=100)
remodel.fit(np.array([x]).T)
remodel.emissionprob_ # ダイスの面の出現確率は揃っていて、不正を見抜くことができない
array([[ 0.165,  0.17 ,  0.17 ,  0.165,  0.165,  0.165]])
# 復数のダイスがあると仮定してみる

remodel = hmm.MultinomialHMM(n_components=2, n_iter=100) # 2個のダイスを想定する
remodel.fit(np.array([x]).T)
remodel.transmat_
array([[ 0.962044  ,  0.037956  ],
       [ 0.18593815,  0.81406185]])
remodel.emissionprob_
array([[  1.31447796e-01,   2.10098809e-01,   1.31451592e-01,
          2.03925627e-01,   1.27562672e-01,   1.95513504e-01],
       [  3.07222812e-01,   2.71009724e-05,   3.33400979e-01,
          4.49521672e-08,   3.23691279e-01,   3.56577840e-02]])
log_p, state_seq = remodel.decode(np.array([x]).T)
state_seq # 綺麗な結果
array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
       0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0])
# ただし、事前知識がないのでfit()の結果は安定しない…
# 次のような結果によくなる
remodel.fit(np.array([x]).T)
log_p, state_seq = remodel.decode(np.array([x]).T)
state_seq
array([0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0,
       1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0,
       0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,
       1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1,
       0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1,
       1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1,
       0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0,
       1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0,
       1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1])
# このようなときは、初期値を明示的に設定してみる

remodel = hmm.MultinomialHMM(n_components=2, n_iter=10000, init_params="st", params="ste")
remodel.emissionprob_ = np.array([[1/6, 1/6, 1/6, 1/6, 1/6, 1/6],
                                  [2/6,   0, 2/6,   0, 2/6,   0]]) # 奇数のみ
remodel.fit(np.array([x]).T)
log_p, state_seq = remodel.decode(np.array([x]).T)
state_seq
/opt/conda/lib/python3.6/site-packages/hmmlearn/hmm.py:405: RuntimeWarning: divide by zero encountered in log
  return np.log(self.emissionprob_)[:, np.concatenate(X)].T





array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

参考:

あやしいサイコロと『隠れマルコフモデル』 | 株式会社フォワードネットワーク

API Reference — hmmlearn 0.2.1 documentation