ユニファ開発者ブログ

ユニファ株式会社システム開発部メンバーによるブログです。

保育士さんの行動推移モデルで理解するマルコフ連鎖

はじめに


こんにちわ、研究開発部の島田です。今回はマルコフ連鎖についてのお話です。

マルコフ連鎖は様々なところで応用されており、イメージしやすい例だとWeb広告のアトリビューション解析などでしょうか。これは、ユーザーがインターネット上で何かを購入するまでの間にどんなWeb広告に接触してどんな行動を取るのか、こういった行動ログをマルコフモデルに当てはめることでどのWeb広告が価値が高いかがわかるようになります。他にもベイズ統計学や強化学習などの分野でもマルコフモデルの考え方が重要になってきます。

名前だけ聞くとすごく難しそうなこのマルコフ連鎖ですが、考え方自体はシンプルで理解しやすいです。 今回はこのマルコフ連鎖について理解し、簡単なシミュレーションも行ってみたいと思います。

マルコフ連鎖とは


マルコフ連鎖についてWiki先生のご説明を見てみます。

マルコフ連鎖(マルコフれんさ、英: Markov chain)とは、確率過程の一種であるマルコフ過程のうち、とりうる状態が離散的(有限または可算)なもの(離散状態マルコフ過程)をいう。また特に、時間が離散的なもの(時刻は添え字で表される)を指すことが多い(他に連続時間マルコフ過程というものもあり、これは時刻が連続である)。マルコフ連鎖は、未来の挙動が現在の値だけで決定され、過去の挙動と無関係である(マルコフ性)。各時刻において起こる状態変化(遷移または推移)に関して、マルコフ連鎖は遷移確率が過去の状態によらず、現在の状態のみによる系列である。特に重要な確率過程として、様々な分野に応用される。

引用元:マルコフ連鎖 - Wikipedia

何言ってるか全くわからんよ・・・Wiki先生。。。

ということで、このブログではまず数学的に小難しいことは無視してイメージだけ掴める様にしてみます。 その後一般的な解釈を数式も交えて理解していきたいと思います。

イメージで理解してみるマルコフ連鎖


マルコフ連鎖をイメージで掴むために、今回は保育士さんの行動に当てはめて理解してみたいと思います。

ある保育士さんが「職員室・園庭・保育室」のいずれかに滞在しているとします。 実際はもっと色々な場所を行ったり来たりして、更に他の要因も複雑に絡み合ってるので簡単な確率モデルでは表せないと思いますが、ここではマルコフ連鎖を理解することが目的なので簡単なモデル(「保育士さん行動推移モデル」と呼ぶ)で考えます。

例えば、今ある保育士さんが保育室に園児たちと一緒にいた所、他の保育士さんからのコールなどで10分後に職員室に行く確率が50%だったとします。職員室での仕事を済ませ、次の10分後はそのまま職員室に滞在する確率を20%、園庭に行く確率を80%とします。そして次に・・・・といった様に複数の状態を時系列的に行き来すること、これこそがマルコフ連鎖です。

f:id:unifa_tech:20190920093726p:plain:w350
(保育士女の子イラスト 出典: かわいいフリー素材集 いらすとや

ここで大事なこととして、マルコフ連鎖は現在の状態だけが次の状態を決めます。先ほどの保育士さんの例で言うと、園庭に行くかどうかはあくまで職員室での状態だけが関係しています。

保育士さん行動推移モデルについてもう少し一般的な解釈


先ほどの説明で何となくマルコフ連鎖についてイメージがつきました。ここからは数式も交えながらもう少し一般的な話をしていきたいと思います。

θを確率変数とすると、下記の様な状態変化の推移のことを確率過程と言います。

θ_{1}→θ_{2}→・・・θ_{t-1}→θ_{t}

この状態変化の推移において、新しい状態θ_{t}が直前の状態θ_{t-1}だけから決定する様な性質をマルコフ性と呼びます。そしてこのマルコフ性を有する確率過程のことをマルコフ連鎖と言います。 一般的には下記の様な式で表せます。

f(θ_{t}|θ_{t-1},・・・,θ_{2},θ_{1}) = f(θ_{t}|θ_{t-1})

では、今回考えている保育士さんの行動推移モデルについて下図のように確率も定義してみます。下図は状態遷移図とも呼ばれ、各場所から場所へ移る確率(遷移核と言う)を示しています。

f:id:unifa_tech:20190920100341p:plain:w400

例えば今ある瞬間に職員室(S)に滞在していたとすると、そこから各場所に移動する確率は下記のようになります。

S → H : 60%, S → E : 30%, S → S : 10%

※S → Sは、職員室に滞在し続けるという意味です。

いちいちこの場所から次の場所への移動確率は○%、と書くのは非常に手間なので遷移核Pは一般的には行列を使って以下のように表します。

$$ P= \begin{pmatrix} 0.1 & 0.6 & 0.3 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.5 & 0.3 \end{pmatrix} $$

保育士さん行動推移モデルから保育士さんの滞在する場所を見積もる


初期状態として、ある保育士さんが職員室にいたとします。その時の確率分布は下記のように書けます。 $$ P(0)= \begin{pmatrix}1&0&0 \end{pmatrix}\ $$

20分後に保育士さんはどこにいる確率が高いかを遷移核行列を使って計算してみます。

$$ P(2)= \begin{pmatrix}1&0&0 \end{pmatrix}\ \begin{pmatrix} 0.1 & 0.6 & 0.3 \\ 0.3 & 0.4 & 0.3 \\ 0.2 & 0.5 & 0.3 \end{pmatrix}^2= \begin{pmatrix}0.25&0.45&0.3\end{pmatrix}\ $$

これは「現在いる場所がわかった時に20分後にどこにいるか」を示しています。結果を見てみると、保育士さんが「現在職員室にいた時に20分後に保育室に行く確率は45%」であることがわかります。

次に時刻変化に対する確率分布がどう変わっていくかを見てみます。 さすがにこれ以上の行列計算を手計算で行うのはシンドイので、コンピュータシミュレーションさせてみます。

今回は統計解析に非常に優れたプログラミング言語であるRを使ってみたいと思います。また、Rにはマルコフ連鎖を便利に使えるパッケージもあるのでそちらも使ってみます。このくらいのシミュレーションであれば実際にはパッケージを使わなくても簡単に出来ちゃいますね。

library(markovchain)

teacherStates <- c("syokuinsitsu", "hoikusitsu", "entei")
byRow <- TRUE
locateMatrix <- matrix(data = c(0.1, 0.6, 0.3,
                               0.3, 0.4, 0.3,
                               0.2, 0.5, 0.3), byrow = byRow, nrow = 3)
mcState <- new("markovchain", states = teacherStates, byrow = byRow,
               transitionMatrix = locateMatrix, name = "teacher_locate")

# set of the initial state is "syokuinsitsu"
initialState <- c(1, 0, 0)

iterlation <- 13
proba_s[1] <- c(initialState[1])
proba_h[1] <- c(initialState[2])
proba_e[1] <- c(initialState[3])
for (i in 1:iterlation)
{
  proba_matrix <- initialState * (mcState ^ i)
  proba_s[i+1] <- proba_matrix[1]
  proba_h[i+1] <- proba_matrix[2]
  proba_e[i+1] <- proba_matrix[3]
}

# make glaph
plot(proba_s, type="o", col="red", xlab="step", ylab="probability", ylim=c(0,1), pch=16)
par(new=T)
plot(proba_h, type="o", xlab="", ylab="", ylim=c(0,1), axes=FALSE, main="", col="blue", pch=17)
par(new=T)
plot(proba_e, type="o", xlab="", ylab="", ylim=c(0,1), axes=FALSE, main="", col="green", pch=15)
legend(x=2,y=1,legend = teacherStates, pch=c(16,17,15), col=c("red","blue","green"))

結果が下の図です。50分経過したあたりから確率分布に変化が無くなっているのがわかります。これを定常分布への収束と言ったりします。つまり今回の設定の場合、職員室に滞在していたある保育士さんが十分に時間が経った時には保育室に滞在している可能性が一番高いことがわかります。

f:id:unifa_tech:20190920135638p:plain

最後に


マルコフ連鎖についてイメージでわかるように簡単な説明とシミュレーションを行なって結果を見てみました。

今回は遷移核がわかっている状態で定常分布を求めましたが、現実問題は定常分布がわかっていて遷移核がわからないケースがほとんどだと思います。そういった場合は、マルコフ連鎖の定常分布に収束するための条件(詳細釣り合い)を考えなくてはなりません。この考えを発展させていくとマルコフ連鎖モンテカルロ法(MCMC)へと繋がっていくわけですね。

最後に、今回は初期状態を職員室に滞在していることとしてシミュレーションしましたが、初期状態を変えたら結果はどうなるのでしょうか?例えば初期状態を保育室にした場合や、職員室50%,保育室50%,園庭0%といった具合に。 もし興味がある方は是非トライしてみてください。5回分の行列積計算であれば頑張れば手計算でも出来ると思います!

参考