黄昏より暗きもの、血の流れより赤きもの

黄昏より暗きもの、血の流れより赤きもの

自分の好きな事を好きなように書いて行きます。

【Advent Calendar 12月12日版 】Mathematica(Raspberry Pi版)のMatrixPower関数を使い、行列のn乗(累乗)を計算してみた【日曜数学】

ご挨拶

Advent Calendar 12/12は私@web_studying(ryuichi69)が書かせて頂きます。Raspberry Piというコンピューターに搭載されている計算処理ソフトMathematicaを使った記事を書こうと思います。Mathematica(Raspberry Pi版)と言う事で、「Mathematica」「Raspberry Pi」の双方のadvent calendarに登録させて頂きました。

さらにAdventerの「日曜数学 Advent Calendar 2015」の方にも出典。id:tsujimottor氏のブログを良く読ませて頂いているのもあり、欲張って合計3つコーナーに出典する結果となりました。

今日は自分が持っているマイコン「Raspberry Pi」にバンドルされているMathematicaの話をしようと思います。Raspberry Piと言えば、ブレッドボードに差し込んでLEDを点灯させたり、リモコンの発信器をつないでTVリモコンの代わりにしたりと言う動きが盛んですね。今日はそんなRaspberry Piの使い方の一例を紹介致します。

Raspberry PiはMineCraftのようなゲームも入っていますし、大学の研究機関向けの数学計算ソフト「Mathematica」がバンドルされています。今日はこれの使い方をと思い、大学受験で言うところの確率漸化式、又は大学の確率論の所で勉強するマルコフ過程という分野を通して、行列のn乗を求めていきます。

マルコフ過程の応用例

マルコフ過程は、ある数量を予測したり、ある情報の重み付けをするための手段として用いられます。例えば1年前の天気や人口などの情報から、一つ前の状態一つ前の状態から現在の状態への遷移に注目し、そこから現在の情報を推測したり、重要な情報かを調べたりするときに使われます。

Webページへの応用例は、googlePageRank[1][2]などがマルコフ過程を用いてページのランク付けを行なった例があります。以下本文を読むのに必要な用語を整理していきます。

本文を読むのに必要な確率用語

確率変数

例えば1〜6が書かれたカードが1枚ずつあり、2回引く事を考えます。その2回の出目をそれぞれX(1),X(2)とします。このとき出目X(1)、X(2)のように試行Tの結果として得られる変数を確率変数と言います。

条件付き確率

このとき1回目に1の出目が出る確率は、以下のように表記します。


f:id:program_study:20151122165821p:plain

さらに1回目に1のカードを引いた条件の元で、2回目に2カードを引く確率と言うように、ある条件の下で起こりうる確率の事を条件付き確率と言います。これを記号を使って表すと以下のようになります。


f:id:program_study:20151122165842p:plain

本日話すマルコフ過程の場合ですと、時刻nの時の状態1から、時刻(n+1)に状態2に遷移する条件付き確率P_{12}=P(X_{n+1}=2|X_{n}=1)を求めるときなどに関連します。

確率過程とマルコフ過程

次に確率過程とマルコフ過程について説明します。nを0以上の整数とします。ある時刻nに応じ、確率変数が変化する過程の事を確率過程と言います。このとき時刻(n+1)における確率変数が、時刻nの状態に依って決まる確率過程の事をマルコフ過程と言います。

例題

nを0以上の整数とし、下図のような三角形ABCと粒子Pがある。時刻n=0にPは頂点Aにいるものとし、AとB、BとC、CとAの間を互いに0.5の確率で動きます。時刻n=10のときにPが頂点A、B、Cに居る確率をそれぞれ求めてください。


f:id:program_study:20151122165344j:plain

解答

Mathematicaに至るまで

まず時刻nにPがA、B、Cに居る確率をそれぞれa_nb_nc_nとおきます。まず時刻nのときにPがAにいて、時刻(n+1)にPがBに移る条件付き確率P_{AB} = 0.5となります。同様にP_{BC} = 0.5P_{CA} = 0.5となります。さて、時刻nのPの位置のみに着目し、時刻(n+1)のPの位置についての漸化式を立てると、以下の式1のようになります。計算のため、時刻n=0の時の状況を式2のようにまとめておきます。


f:id:program_study:20151122165456p:plainf:id:program_study:20151122165530p:plain


上記を行列とベクトルに直すと式3が得られます。


f:id:program_study:20151122165554p:plain


ここで上記の3行3列の行列のように、次の状態に推移することを表した行列の事を推移確率行列と言います。この推移確率行列をPとおき、他の二つのベクトルをA_{n+1},A_nとおきます。このとき以下の式4のようになります。


f:id:program_study:20151122165633p:plain

上記の式から分かる通り、初期状態と行列のn乗を求める事で、時刻nのときにPがA,B,Cのいずれかに居る確率を求めることができます。(この事情は一般化され、チャップマン=コルモゴロフ方程式[3]という定理にもなっています。詳細を知りたい方は参考文献の[3]を参照して下さい。)

なお今回はn=10の状況を考えれば良く、以下の式6のようになります。


f:id:program_study:20151122165700p:plain

Mathematicaで計算

依ってこの行列Pの10乗を求めてやります。尚計算する行列が正方行列の場合に限り、MatrixPower関数を使ってやります。例えば今回のように3行3列の行列の10乗を求める場合、以下のように入力します。

In[10]:= MatrixPower[{{0,0.5,0.5},{0.5,0,0.5},{0.5,0.5,0}},10]
Out[10]= {{0.333984,0.333008,0.333008},{0.333008,0.333984,0.333008},{0.333008,0.333008,0.333984}}

ここでMatrixPower[{正方行列},次数]の順に引数を入力します。第一引数の正方行列を1行ずつカッコでかこって入力して行く事にご注意下さい。ちなみに「.」を使うと行列のかけ算もできるようです。

In[11]:= {MatrixPower[{{0,0.5,0.5},{0.5,0,0.5},{0.5,0.5,0}},10]}.{1,0,0}
Out[11]= {{0.333984,0.333008,0.333008}}

ここでa[10]=0.333984,b[10]=0.333008,c[10]=0.333008となるので、時刻n=10にPがA、B、Cに居る確率はそれぞれ0.333984,0.333008,0.333008となります。

おまけ:桁落ちを確かめる

又小数点の桁落ちが無いか確かめます。まず時刻nにおいて点PはA,B,Cのいずれかの位置に居ます。この事から式を立ててやると、


f:id:program_study:20151122165730p:plain

となります。これはa_nとb_nとc_nの総和が1と言う事を表しています。n=10なので先程求めた答えをa[10]+b[10]+c[10] = 0.333984 + 0.333008 + 0.333008 = 1.000000。これをみるに行列で計算しているにもかかわらず桁落ちしてないのも、mathematica凄いな〜って思います(数値計算のプログラム組む時って、この辺が難しいですからね)

是非合わせて読んでみて下さい!

[1]PageRankの基本概念:http://tnt.math.se.tmu.ac.jp/labo/grad/2004/masa/graph/6.html
[2]PageRankアルゴリズムおよびそれに関連する研究について 京都大学 総合人間学部 認知情報学系 数理情報論分野 宮嶋健人:http://www.kentmiyajima.com/document/pagerank.pdf
[3]Aritalab:Lecture/NetworkBiology/Markov Chains:http://www.metabolomics.jp/wiki/Aritalab:Lecture/NetworkBiology/Markov_Chains

後書き

本記事ではMathematicaの特徴を、実際の数学の問題を解きながら説明してみました。手計算の結果a[n]が1/3に収束することが分かっていたので、0.3333…と3だけが出ないような時刻nに設定しました。問題作るのってバランス調整が難しいです。

そして本問題は行列以外に簡単な漸化式などで解く事もできますので、是非とも違う関数を作ってMathematicaの練習をしてみて下さい。最後になりますがadvent calenderを書いている皆さん、記事の方がんばって下さい。