第 4 章 状態空間モデルによる分析 21
4.3 ベイズ法によるパラメータ推定
ここまでで実験データを取得し、データから分析方針を立て、状態空間モデルを構築し た。構築した状態空間モデルには、𝑘1 や𝑘2などの推定しなければならないパラメータが ある。パラメータ推定の際、尤度関数が複雑であることやパラメータ数の多さの問題など から、最尤推定では困難が生じる。そこで、パラメータを推定するために、ベイズ統計学 を基盤とする推定方法であるベイズ法を用いる。ベイズ法による分析やパラメータ推定は、
次のような手順で進める[9]。
1. データやその収集過程に関する知識に基づき、確率モデルを作る
2. 観測されたデータを条件として、モデルのパラメータ事後確率分布を計算する 3. データに対するモデルの当てはまりや、計算された事後確率分布の評価をおこなう 手順1については前節で終えたので、以降では残りの手順2と手順3を説明する。ここ
では、Stan[11]というソフトウェアを使用する。Stanは手順2と手順3の計算を自動化す
るために用いられる。本節では、はじめに手順2と手順3の、パラメータの推定方法とそ の評価方法について簡単に説明する。つづいて、Stanを用いて計算を実行する際の準備に ついて説明する。
4.3.1
パラメータの事後分布
ベイズ統計学では、モデルのパラメータは何らかの確率分布に従う確率変数であるとす る。データが得られる前に想定したパラメータの確率分布を事前確率分布といい、実際に 得られたデータで条件づけられた確率分布を事後確率分布という。事前確率分布と事後確 率分布を、単に事前分布と事後分布と呼ぶことも多い。ベイズの定理より、事後分布は、
事前分布とデータから計算される尤度の積に比例する。すなわち、事前分布に想定する分 布によって、事後分布に反映される尤度の影響の大きさが変わってくるのである。対象に
対する前提知識や想定がない場合、事前分布に無情報事前分布を用いることが多い。無情 報事前分布とは、パラメータが取り得るすべての値に対して、みな等しく平等に低い確率 密度を割り当てた確率分布である。たとえば、裾の広い一様分布や分散が非常に大きい正 規分布などが無情報事前分布として用いられることが多い。
このように、パラメータが事前分布と尤度の積に比例する事後分布に従うという想定は、
次の3つのようなことが利点として挙げられる。1つ目は、事後分布を推定した後に新た なデータが入手できた時、推定した事後分布を新たな事前分布として、再び事後分布を計 算することができることである。これをベイズ更新と呼ぶ。2つ目は、現象の観察や理論 的背景から得られる想定を事前分布に反映させることができることである。3つ目は、パ ラメータが取りうる値の範囲を確率的に知ることができることである。
本研究で作成したモデルのパラメータは、各粒径域𝑖（最大粒径域を除く）におけるミ カエリス・メンテン式のパラメータ𝑘1 と𝑘2、放流率の𝑞、および各粒径域の割合の初期 値x1である。初期値x1は、反応槽の真の粒径域ごとの割合の初期値であり観測すること ができないので、パラメータとして推定する。構築したモデルでは粒径域の数を11とし ているため成長割合𝑝𝑖 は10個であり、パラメータ数は全部で32個である。また、事前 分布には無情報事前分布として一様分布を使用するが、パラメータの取りうる値には上限 と下限に対して制限を設けている。詳しくはStanの準備にて後述する。
4.3.2 MCMC
法による計算
MCMCとはMarkov Chain Monte Carlo（マルコフ連鎖モンテカルロ）の略称で、これ は乱数発生のためのアルゴリズムである。この乱数発生アルゴリズムであるMCMCを事 後分布の計算に応用したものを、ここではMCMC法と呼ぶことにする。なぜ事後分布の 計算にMCMCを応用するかは、次のような理由からである。ベイズの定理より、事後分 布の確率密度関数は、尤度関数と事前分布の確率密度関数の積を、尤度関数と事前分布の 確率密度関数の積を積分したもので除したものである。しかし、この事後分布の確率密度 関数は、モデルが複雑となることが多く、特にパラメータが多い時には非常に複雑になり、
積分計算が困難になる。そこで、MCMCによって事後分布の確率密度関数に従う乱数を 発生させることで、パラメータがある範囲の値を取る確率や、事後分布の期待値などが計 算できるようになる。
MCMC法による計算方法には、MH法（メトロポリス・ヘイスティング法）やHMC法
（ハミルトニアン・モンテカルロ法）などがある。MH法とは「仮の数値をそれぞれの変数 に投入し、他の変数値を固定したときにある変数の2つの数値を比較して、データの当て はまり（尤度、もっともらしさを表す統計値）と事前分布の確率の積が高い方を採用して いく手順を取って、変数値の分布（事後分布）を描き出す方法である（計算の過程で、あ る一定の割合で尤度と事前確率の積が低い値も採用するようになっている）」[12]。HMC 法の基本的な考え方はMH法と同じであり、Stanで使用するのはHMC法である。
MCMC法による計算の結果、抽出される値がある事後分布に収束しているかどうかを 判断するために、収束診断が必要である。この収束診断法方法は複数あるが、本研究にお いては次の2つの方法で収束しているかどうかの判断をおこなった。1つは、トレースプ ロットと呼ばれる時系列プロットの目視による判定である。MCMC法では初期値の異な る複数の連鎖（chains）で乱数を発生させることができ、これを時系列でプロットしたも のがトレースプロットである。トレースプロットですべての連鎖が別れることなく交わっ ていることが収束の判断基準である。もう1つはRhatが1.1未満になっているかどうか である。Rhatは、事後分布からのサンプリングが収束したかを見る指標で、1へ収束する 指標である。初期値の異なる連鎖の収束先が異なっていた場合などは、Rhatは大きくな る。すべての𝑘1と𝑘2、𝑞、および初期値x1について、Rhatが1.1未満であることが収束 の判断基準である。
4.3.3 Stan
の準備
MCMC法を実行するためのソフトウェア環境にはBUGSやStanなどがあるが、本研 究ではStanを利用する。Stanに受け渡すためのデータは、すべてRで準備している。R で生データが保存されているサーバーのデータベースから必要なデータを抽出し、モデル の推定に必要なデータセットを作成した。StanはRStan[13]というパッケージを用いて、
Rから呼び出している。
MCMC法の計算にあたり、パラメータそれぞれに次のような制約を設けた。放流率𝑞 に関しては、排液によって反応槽から結晶がなくなる割合なので下限を0とし、上限を 0.4とした。上限を0.4としたのは、単位時間あたりに、反応槽にある結晶がほぼ半分な くなることはないと考えられるためである。ミカエリス・メンテン式のパラメータ𝑘1 は、
溶液濃度なので下限を0とし、上限は1とした。3章で述べたように、本実験装置におい
ては0.004mol/kg程度が実験に適した濃度の上限であるため、上限を1としても推定に支
障はないと考えられる。ミカエリス・メンテン式のパラメータ𝑘2は、成長割合（の上限）
なので下限を0とし、上限は1とした。
つぎに、MCMC法を実行する際の計算回数などの指定について説明する。MCMC法で は、アルゴリズムの都合上、乱数が初期値へ依存してしまうことや、連続する採用値に自 己相関が含まれてしまうことが起こりえる。こういった問題を回避するために、次のよう な対策をする。はじめにウォームアップ期間（warmup）を設定する。ウォームアップ期 間とは、生成された乱数を切り捨てる期間である。ウォームアップ期間を指定することで、
初期値への依存性を下げることができる。つぎに間引き間隔（thin）の設定をする。間引 き間隔とは、生成された乱数を間引く間隔の指定である。たとえば間引き区間を2とする と、生成された乱数のうち1つ目の値を採用し、2つ目の値を切り捨て、3つ目の値を採 用するといったようになる。乱数を間引いて採用することで、乱数の自己相関が起きるこ とを防ぐ。さいごにマルコフ連鎖の数（chains）を設定する。マルコフ連鎖の数とは、初 期値の異なる乱数の連鎖を発生させる数の指定である。異なる初期値から出発した乱数の 連鎖が交わっているかどうかを確認することによって、事後分布が収束したかどうかの目 安となる。
本研究におけるMCMC法の計算は、20万回おこない、初期値からの10万回を切り捨 て（warmup=100,000）、以後の100回ごとに結果を抽出（thin=100）する計算を初期値の 異なる3本の連鎖（chains=3）についておこない、事後分布を求めた。すなわち結果とし て得られる事後分布からの標本の大きさは3千である。なお、サンプリング結果が再現で きるよう乱数のシードを1229とした。