orizuru

つながる.見える.わかる IoTソリュ-ション

階層ベイズモデル

約 13 分
階層ベイズモデル

はじめに

 書籍「ベイズモデリングの世界」に「階層モデルで『個性』をとらえる」という解説がある。著者は北大の生態学を専門とする久保拓弥氏である。架空の植物を考え、その結実確率(胚珠が種子になる確率)を階層ベイズモデルを用いて推定する過程が分かりやすく解説されている。階層ベイズモデルのパラメータを推定する代表的な手法には、経験ベイズ法とMCMC法がある。久保氏の解説では前者を用いた推定がRを用いて行われている。MCMC法による解法の詳細には言及していないが、Rのサンプルプログラムへのリンクが記載されている。今回は、このサンプルプログラムと同じことを、前回紹介したPyMCを用いてなぞってみたい。今回使うのはPyMC3である(前回はPyMC2であった)。

ソースコード

 今回実装したコードはここにあるsample_0.ipynbである。

問題設定

 ここで取り上げる問題は以下の通り。

  • 架空植物のある1個体を選択したとき、その植物がいくつ種子を付けるのかを推定したい。
  • この植物は胚珠(種子の元になる器官)を必ず10個持つ。
  • 全ての胚珠が必ず種子になるわけではない。ある個体の種子は4個、別の個体は0個だったりする。
  • この架空植物を100個体観測する。この観測データから各個体の結実確率を推定したい。

  • 全ての個体が10個の種子を作れば、全種子数は1000(=10\times100)、全ての個体がどれも種子を作らなければ全種子数は0(=0\times100)になる。データはここからダウンロードできる。データの中身は以下の通り。

    “plant.ID”,”y”,”alpha”
    1,0,-4.15958227809385
    2,0,-3.71836219544951
    3,0,-3.06901965190728
    4,0,-3.05421060686106
    5,0,-2.88065618323557
    6,2,-2.70559527531842
    7,1,-2.44170265669845
    8,1,-2.22132654301567

    解析に使う列はyである。この列には各個体の種子数が記載されている。横軸に種子数y、縦軸に個体数pを取ると以下のグラフを得る。

    種子数0の個体数は7、種子数1の個体数は10であることなどが分かる。種子数0から10までの個体数を全て足すと100になることに注意する。

    モデル化

     ある1個体を考える。この個体が種子をy個作る確率をモデル化したい。各個体には胚珠(種子の元になる器官)が10個あり、各胚珠は種子になるかならないかのどちらかである。従って、次の二項分布{\rm Bin}で表すことができる。

    (1)    \begin{eqnarray*} p(y|q)&=&{\rm Bin}(y|10,q) \\ {\rm Bin}(y|10,q)&=&{}_{10} C _y\;q^y(1-q)^{10-y} \end{eqnarray*}

    ここで、qは胚珠が種子になる確率(結実確率)である。このとき、全ての観測データが実現する確率は

    (2)    \begin{equation*} p(y_1,\cdots,y_{100}|q_1,\cdots,q_{100})=\prod_{i=1}^{100}p(y_i|q_i) \end{equation*}

    となる。

    個体差を考慮しない予測

     いま結実確率q_iが個体に寄らず一定であると仮定すると、対数尤度は

    (3)    \begin{eqnarray*} \ln{p(y_1,\cdots,y_{100}|q) &=& \sum_{i=1}^{100}\ln{p(y_i|q)}\\ &=& \sum_{i=1}^{100}\ln{{\rm Bin}(y_i|10,q)}\\ &=& \sum_{i=1}^{100}\left[y_i\ln{q}+(10-y_i)\ln{(1-q)} \right] + {\rm const.} \end{eqnarray*}

    となる。この式をqで微分することにより最尤推定値\hat{q}を得る。

    (4)    \begin{equation*} \hat{q}=\frac{1}{1000}\sum_{i=1}^{100}y_i \end{equation*}

    観測値を用いて\hat{q}を計算すると0.496となる。このときの種子数と個体数の間の関係式は

    (5)    \begin{equation*} p=100\;{\rm Bin}(y|10,\hat{q}) \end{equation*}

    となり、対応するグラフは以下の点線となる。

    明らかに観測データを再現できていない。

     ここまでの解析をPyMC3で行う。コードは以下の通り。

  • 3行目:qに対して事前確率を定義する。lower=0からupper=1の範囲で定義される一様分布とした。
  • 4行目:先に示した尤度(二項分布)を定義する。observedに観測値をysを与える。
  • 5行目:事後確率最大化(MAP推定)を行う。

  • qymodelが必要とするstochastic変数である。PyMC2では、ユーザが明示的にこれらをModelのコンストラクタに渡していた。PyMC3では、上のコードのように、withの直下にモデルに必要な各種変数を記述すれば良い。上記の出力は以下の通り。

    0.4960000053023121

    先に示した値と同じである(qの事前分布として一様分布を与えたから、MAP推定値と最尤推定値は一致する)。

    個体差を考慮する予測

     次に結実確率qを個体iごとに異なるものとする。

    (6)    \begin{eqnarray*} p(y_i|q_i) &=& {\rm Bin}(y_i|10,q_i) \\ &=& {\rm Bin}(y_i|10,f_S(\beta+\alpha_i))\equiv p(y_i|\beta,\alpha_i) \\ f_S(z)&=&\frac{1}{1+\exp{(-z)}} \end{eqnarray*}

    ここで、f_S(z)はシグモイド関数である。\betaは個体に依存しないパラメータを、\alpha_iは個体に依存するパラメータを表す。\beta\alpha_iについて次の事前確率を定義する。

    (7)    \begin{eqnarray*} p(\beta|\tau_{\beta})&=&\mathcal{N}(\beta|0,1/\tau_{\beta}) \\ p(\alpha_i|\tau_{\alpha})&=&\mathcal{N}(\alpha_i|0,1/\tau_{\alpha}) \end{eqnarray*}

    ここで、\mathcal{N}(x|\mu,1/\tau)は平均\mu、精度\tauの正規分布を表す(分散\sigma^2の逆数を精度(precision)と呼ぶ)。個体に依存しない\betaについては分散の大きな(1/\tau_{\beta}=100)分布を考える。一方、個体に依存する\alpha_iに対しては、もう少し注意深く考察したい。そのため、次の事前確率を定義する。

    (8)    \begin{eqnarray*} p(\tau_{\alpha}|a,b) &=& {\rm Gam}(\tau_{\alpha}|a,b) \\ {\rm Gam}(x|a,b) &=& \frac{b^a}{\Gamma(a)}x^{a-1}\exp{(-bx)} \end{eqnarray*}

    ここで、{\rm Gam}(x|a,b)はガンマ分布であり、これは正の実数を生成する。\tau_{\alpha}を定数ではなく確率から生成される値とすることで、観測データに適応させて\alpha_iを決めることができる(ただし、abはあらかじめ決めておくハイパーパラメータである)。このように事前分布のパラメータに対し、さらに事前分布を導入するモデルを階層ベイズモデルと呼ぶ。\vec{y}=(y_1,\cdots,y_{100})\vec{\alpha}=(\alpha_1,\cdots,\alpha_{100})と置くと、求めたい事後分布は次式で与えられる。

    (9)    \begin{equation*} p(\vec{\alpha},\beta,\tau_{\alpha}|\vec{y}) = \frac{ \prod_{i=1}^{100}\;p(y_i|\beta, \alpha_i)p(\beta|\tau_{\beta})p(\alpha_i|\tau_{\alpha})p(\tau_{\alpha}|a,b)} {p(\vec{y})} \end{equation*}

    ここで、右辺分母は次式で計算される。

    (10)    \begin{equation*} p(\vec{y})=\prod_{i=1}^{100}\int d\beta\;d\alpha_i\;d\tau_{\alpha}\;p(y_i|\beta, \alpha_i)p(\beta|\tau_{\beta})p(\alpha_i|\tau_{\alpha})p(\tau_{\alpha}|a,b) \end{equation*}

    この積分は解析的に計算できないので、ここから先の計算はPyMC3を用いてMCMC法で近似する。コードは以下の通り。

  • 3行目:p(\beta|\tau_{\beta})を定義する。tau=\tau_{\beta}=0.01
  • 4行目:p(\tau_{\alpha}|a,b)を定義する。alpha=a=0.01,beta=b=0.01
  • 5行目:p(\alpha_i|\tau_{\alpha})を定義する。shape=100として、100個の\alpha_iを生成する。
  • 6行目:f_S(\beta + \alpha_i)を計算する。
  • 7行目:p(y_i|\beta,\alpha_i)を定義する。observedには観測値ysを与える。
  • 8行目:MAP推定を行う。
  • 9,10行目:サンプリングを行う。

  • 上のコードを実行したあと、以下のコードで収束しているか否かを判定する。

    出力は以下の通り。

    {‘alpha’: array([0.99992116, 0.99997074, 0.99997656, 1.00101515, 1.0000671 ,
    0.99997194, 0.99990008, 1.00046983, 0.99990016, 1.00004488,
    0.99997887, 0.99996494, 1.00001114, 0.99993037, 1.00009044,
    0.9999153 , 0.99990012, 0.99999316, 0.99990019, 0.99990849,
    0.9999183 , 1.00019395, 0.99997166, 0.99990267, 1.00025804,
    1.00006084, 1.00003207, 1.00001905, 1.00002649, 0.99999359,
    0.99991317, 0.99990444, 1.00010107, 1.00044065, 1.00008226,
    0.9999002 , 0.99991316, 0.99995332, 1.00016158, 1.00016238,
    0.99991051, 0.99994176, 1.00003935, 0.99997246, 1.00052741,
    1.00017737, 0.99995378, 1.00002957, 1.00018027, 0.99999548,
    1.00012819, 0.99990123, 0.99993936, 0.99992739, 0.99991631,
    0.9999227 , 0.9999047 , 1.00006704, 0.99997697, 0.99992804,
    0.9999103 , 0.99990415, 0.99991652, 0.99991608, 1.00098608,
    0.99990014, 0.99990024, 0.9999914 , 1.00003664, 0.99991117,
    0.99990698, 0.99997845, 0.99991626, 0.99992411, 1.00014929,
    0.9999008 , 0.9999 , 0.99992404, 0.99990093, 0.99990898,
    1.00012367, 1.00011088, 0.99995555, 0.99991072, 0.99992426,
    1.00033826, 1.00004621, 0.99990063, 0.99991521, 1.00011476,
    1.00005838, 0.99990038, 0.99991504, 0.99990287, 0.99991835,
    1.00019647, 0.99993153, 0.99993019, 0.99995296, 1.00004081]),
    ‘beta’: 0.9999865601392494,
    ‘tau’: 1.0006839287826985}

    100個の\alpha_i\beta\tau_{\alpha}のGelman-Rubin統計量\hat{R}が出力される。詳細は省くが、経験則から\hat{R}<1.1のとき収束したとみなして良い。上の結果を見ると全変数が1.1未満となっていることが分かる。ここまでの作業で、式(9)で定義される事後確率が求まったことになる。この事後確率を用いると、予測される種子数y^*の分布は次式で計算される。

    (11)    \begin{eqnarray*} p(y^*|\vec{y}) &=& \int d\beta\;d\vec{\alpha}\;d\tau_{\alpha}\ p(y^*|\beta,\vec{\alpha}) p(\vec{\alpha},\beta,\tau_{\alpha}|\vec{y}) \\ &\simeq&\frac{1}{M} \sum_{(\vec{\alpha},\beta,\tau_{\alpha})\sim p(\vec{\alpha},\beta,\tau_{\alpha}|\vec{y})}p(y^*|\beta,\vec{\alpha}) \end{eqnarray*}

    和は事後確率からサンプリングした点を用いて行えば良い。Mはサンプル点の総数である。この計算を行うコードは以下の通り。

  • 1行目:サンプル点の数を取り出す。
  • 2行目:事後確率を計算する。
  • 4行目:yの事後確率のサンプル点を取り出す。

  • これだけのコードでp(y^*|\vec{y})からサンプリングされた点を得ることができる。 このサンプル点からヒストグラムを作れば、yの予測分布を得ることができる。

    次のコードで描画する。

    以下に結果を示す。

    凡例の意味は以下の通り。

  • 青色の丸:観測値
  • 橙色の丸(点線):結実確率に個体差がないと仮定した場合の分布
  • 緑色の丸:結実確率に個体差があると仮定した場合の分布(平均値)
  • 青色の領域:\sigma(標準偏差)の範囲
  • 黄色の領域:2\sigmaの範囲

  • 結実確率に個体差を考慮することにより、観測データをかなり良く再現できるようになることが分かる。

    まとめ

     今回は、階層ベイズモデルにより、個体差を取り込む手順を示した。尤度p(y|\beta,\alpha_i)に対し、パラメータ\beta\alpha_iの事前確率を仮定する。後者の事前確率に含まれるパラメータ\tau_{\alpha}に対し、さらに事前確率を導入することにより、\alpha_iを単なる主観ではなく観測データに適応させて決めることができるようになる。すなわち、観測データから個体差を抽出することができるのである。
     今回の解析ではPyMC3を用いた。前回使用したPyMC2よりモデルを記述しやすくなっている。また、PyMC3ではバックエンドにTheanoが使われており、速度面でも改善が見られる。ご存知のようにTheanoは最近メンテナンスを停止する旨の告知がなされた。PyMC3の開発チームは、PyMC3のバックエンドとして必要となる機能だけではあるが、メンテナンスを引き継ぐようである。また、それと並行してPyMC4の開発が進められている。こちらのバックエンドはTensorFlow Probabilityなるモジュールを使うようだ。PyMC4のリリースはまだまだ先であり、今後もPyMC3の機能拡張やバグフィックスが続けられるとのことである(引用元)

    About The Author

    IoT/AIソリューション事業部(深層/機械学習・画像処理エンジニア)KumadaSeiya
    深層/機械学習と画像処理などを担当。物性理論で博士号を取得。
    http://seiya-kumada.blogspot.jp/
    https://twitter.com/seiya_kumada

    Leave A Reply

    *
    *
    * (公開されません)