階層ベイズモデル

はじめに

書籍「ベイズモデリングの世界」に「階層モデルで『個性』をとらえる」という解説がある。著者は北大の生態学を専門とする久保拓弥氏である。架空の植物を考え、その結実確率（胚珠が種子になる確率）を階層ベイズモデルを用いて推定する過程が分かりやすく解説されている。階層ベイズモデルのパラメータを推定する代表的な手法には、経験ベイズ法と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で行う。コードは以下の通り。

model = pm.Model()
with model:
    q = pm.Uniform('q', lower=0, upper=1)
    y = pm.Binomial('y', n=10, p=q, observed=ys)
    map_estimate = pm.find_MAP(model=model)
print(map_estimate['q'])

model = pm.Model()

with model:

q = pm.Uniform('q', lower=0, upper=1)

y = pm.Binomial('y', n=10, p=q, observed=ys)

map_estimate = pm.find_MAP(model=model)

print(map_estimate['q'])

3行目： $q$

に対して事前確率を定義する。lower=0からupper=1の範囲で定義される一様分布とした。

4行目：先に示した尤度（二項分布）を定義する。observedに観測値をysを与える。

5行目：事後確率最大化（MAP推定）を行う。

qとyはmodelが必要とする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$ を決めることができる（ただし、 $a$ と $b$ はあらかじめ決めておくハイパーパラメータである）。このように事前分布のパラメータに対し、さらに事前分布を導入するモデルを階層ベイズモデルと呼ぶ。 $\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法で近似する。コードは以下の通り。

model = pm.Model()
with model:
    beta = pm.Normal('beta', mu=0, tau=1.0e-2)
    tau = pm.Gamma('tau', alpha=1.0e-02, beta=1.0e-02)
    alpha = pm.Normal('alpha', mu=0, tau=tau , shape=len(ys))
    ymu = pm.math.sigmoid(beta + alpha)
    y = pm.Binomial('y', n=10, p=ymu, observed=ys)
    start = pm.find_MAP()
    step = pm.NUTS()
    trace = pm.sample(5000, start=start, step=step)

model = pm.Model()

with model:

beta = pm.Normal('beta', mu=0, tau=1.0e-2)

tau = pm.Gamma('tau', alpha=1.0e-02, beta=1.0e-02)

alpha = pm.Normal('alpha', mu=0, tau=tau , shape=len(ys))

ymu = pm.math.sigmoid(beta + alpha)

y = pm.Binomial('y', n=10, p=ymu, observed=ys)

start = pm.find_MAP()

step = pm.NUTS()

trace = pm.sample(5000, start=start, step=step)

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行目：サンプリングを行う。

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

# 収束をチェックする。
pm.gelman_rubin(trace)

1 2	# 収束をチェックする。 pm.gelman_rubin(trace)

出力は以下の通り。

{‘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$ はサンプル点の総数である。この計算を行うコードは以下の通り。

sample_size, _ = trace['alpha'].shape
posterior = pm.sample_ppc(trace, sample_size, model=model)

y_posterior = posterior['y']

sample_size, _ = trace['alpha'].shape

posterior = pm.sample_ppc(trace, sample_size, model=model)

y_posterior = posterior['y']

1行目：サンプル点の数を取り出す。

2行目：事後確率を計算する。

4行目： $y$

の事後確率のサンプル点を取り出す。

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

# sample_sizeの数だけヒストグラムを作る。
pred_ys = np.empty((sample_size, bins))
for i, k in enumerate(y_posterior):
    pred_ys[i] = np.histogram(k, bins)[0]

# 平均をとる。
pred_ys = pred_ys.mean(axis=0)
# 標準偏差を計算する。
pred_std = pred_ys.std(axis=0)

# sample_sizeの数だけヒストグラムを作る。

pred_ys = np.empty((sample_size, bins))

for i, k in enumerate(y_posterior):

pred_ys[i] = np.histogram(k, bins)[0]

# 平均をとる。

pred_ys = pred_ys.mean(axis=0)

# 標準偏差を計算する。

pred_std = pred_ys.std(axis=0)

次のコードで描画する。

#  描画する。
plt.figure(figsize=(10,6))
plt.rcParams['font.size'] = 16
plt.scatter(xs, histo_ys, s=100, label='observed data')

plt.plot(xs, binom_ys, ls='--')
plt.scatter(xs, binom_ys, s=100, label='fixed $q

以下に結果を示す。
<img src="/wp-content/uploads/2018/09/fig4.png" alt="" width="623" height="386" class="aligncenter size-full wp-image-5975">
凡例の意味は以下の通り。
<nl></nl>

 	<li>青色の丸：観測値</li>


 	<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>


 	<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>


 	<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>


 	<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>



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

<h2 class="category_header left"><span class="inner left">まとめ</span></h2>


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

plt.xlabel('$y

以下に結果を示す。
<img src="/wp-content/uploads/2018/09/fig4.png" alt="" width="623" height="386" class="aligncenter size-full wp-image-5975">
凡例の意味は以下の通り。
<nl></nl>

 	<li>青色の丸：観測値</li>


 	<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>


 	<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>


 	<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>


 	<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>



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

<h2 class="category_header left"><span class="inner left">まとめ</span></h2>


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

以下に結果を示す。
<img src="/wp-content/uploads/2018/09/fig4.png" alt="" width="623" height="386" class="aligncenter size-full wp-image-5975">
凡例の意味は以下の通り。
<nl></nl>

 	<li>青色の丸：観測値</li>


 	<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>


 	<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>


 	<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>


 	<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>



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

<h2 class="category_header left"><span class="inner left">まとめ</span></h2>


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

plt.plot(xs, pred_mean, ls='-.', c='g')
plt.scatter(xs, pred_mean, s=100, label='non-fixed $q$(mean)')

plt.fill_between(xs, pred_mean - 2*pred_std, pred_mean + 2*pred_std, color='yellow', alpha=0.2, label='$\pm2\sigma

以下に結果を示す。
<img src="/wp-content/uploads/2018/09/fig4.png" alt="" width="623" height="386" class="aligncenter size-full wp-image-5975">
凡例の意味は以下の通り。
<nl></nl>

 	<li>青色の丸：観測値</li>


 	<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>


 	<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>


 	<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>


 	<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>



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

<h2 class="category_header left"><span class="inner left">まとめ</span></h2>


　今回は、階層ベイズモデルにより、個体差を取り込む手順を示した。尤度[latex]p(y|\beta,\alpha_i)[/latex]に対し、パラメータ[latex]\beta[/latex]と[latex]\alpha_i[/latex]の事前確率を仮定する。後者の事前確率に含まれるパラメータ[latex]\tau_{\alpha}[/latex]に対し、さらに事前確率を導入することにより、[latex]\alpha_i[/latex]を単なる主観ではなく観測データに適応させて決めることができるようになる。すなわち、観測データから個体差を抽出することができるのである。
　今回の解析ではPyMC3を用いた。前回使用したPyMC2よりモデルを記述しやすくなっている。また、PyMC3ではバックエンドにTheanoが使われており、速度面でも改善が見られる。ご存知のようにTheanoは最近メンテナンスを停止する旨の告知がなされた。PyMC3の開発チームは、PyMC3のバックエンドとして必要となる機能だけではあるが、メンテナンスを引き継ぐようである。また、それと並行してPyMC4の開発が進められている。こちらのバックエンドはTensorFlow Probabilityなるモジュールを使うようだ。PyMC4のリリースはまだまだ先であり、今後もPyMC3の機能拡張やバグフィックスが続けられるとのことである<a href="https://medium.com/@pymc_devs/theano-tensorflow-and-the-future-of-pymc-6c9987bb19d5" target="_blank" rel="nofollow noopener noreferrer">（引用元）</a>。
)
plt.fill_between(xs, pred_mean - pred_std, pred_mean + pred_std, color='blue', alpha=0.2, label='$\pm\sigma

以下に結果を示す。
<img src="/wp-content/uploads/2018/09/fig4.png" alt="" width="623" height="386" class="aligncenter size-full wp-image-5975">
凡例の意味は以下の通り。
<nl></nl>

 	<li>青色の丸：観測値</li>


 	<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>


 	<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>


 	<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>


 	<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>



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

<h2 class="category_header left"><span class="inner left">まとめ</span></h2>


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

plt.gca().xaxis.set_minor_locator(tick.MultipleLocator(1))
plt.grid(which='minor')
plt.legend(loc='upper left')
plt.show()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

# 描画する。

plt.figure(figsize=(10,6))

plt.rcParams['font.size'] = 16

plt.scatter(xs, histo_ys, s=100, label='observed data')

plt.plot(xs, binom_ys, ls='--')

plt.scatter(xs, binom_ys, s=100, label='fixed $q

以下に結果を示す。

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

<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>

<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>

<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>

<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>

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

　今回は、階層ベイズモデルにより、個体差を取り込む手順を示した。尤度[latex]p(y|\beta,\alpha_i)[/latex]に対し、パラメータ[latex]\beta[/latex]と[latex]\alpha_i[/latex]の事前確率を仮定する。後者の事前確率に含まれるパラメータ[latex]\tau_{\alpha}[/latex]に対し、さらに事前確率を導入することにより、[latex]\alpha_i[/latex]を単なる主観ではなく観測データに適応させて決めることができるようになる。すなわち、観測データから個体差を抽出することができるのである。

　今回の解析ではPyMC3を用いた。前回使用したPyMC2よりモデルを記述しやすくなっている。また、PyMC3ではバックエンドにTheanoが使われており、速度面でも改善が見られる。ご存知のようにTheanoは最近メンテナンスを停止する旨の告知がなされた。PyMC3の開発チームは、PyMC3のバックエンドとして必要となる機能だけではあるが、メンテナンスを引き継ぐようである。また、それと並行してPyMC4の開発が進められている。こちらのバックエンドはTensorFlow Probabilityなるモジュールを使うようだ。PyMC4のリリースはまだまだ先であり、今後もPyMC3の機能拡張やバグフィックスが続けられるとのことである<a href="https://medium.com/@pymc_devs/theano-tensorflow-and-the-future-of-pymc-6c9987bb19d5" target="_blank" rel="nofollow noopener noreferrer">（引用元）</a>。

)

plt.xlabel('$y

以下に結果を示す。

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

<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>

<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>

<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>

<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>

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

)

plt.ylabel('$p

以下に結果を示す。

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

<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>

<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>

<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>

<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>

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

)

plt.plot(xs, pred_mean, ls='-.', c='g')

plt.scatter(xs, pred_mean, s=100, label='non-fixed $q$(mean)')

plt.fill_between(xs, pred_mean - 2*pred_std, pred_mean + 2*pred_std, color='yellow', alpha=0.2, label='$\pm2\sigma

以下に結果を示す。

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

<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>

<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>

<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>

<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>

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

)

plt.fill_between(xs, pred_mean - pred_std, pred_mean + pred_std, color='blue', alpha=0.2, label='$\pm\sigma

以下に結果を示す。

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

<li>橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布</li>

<li>緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）</li>

<li>青色の領域：[latex]\sigma[/latex]（標準偏差）の範囲</li>

<li>黄色の領域：[latex]2\sigma[/latex]の範囲</li>

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

)

plt.gca().xaxis.set_minor_locator(tick.MultipleLocator(1))

plt.grid(which='minor')

plt.legend(loc='upper left')

plt.show()

以下に結果を示す。

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

青色の丸：観測値

橙色の丸（点線）：結実確率に個体差がないと仮定した場合の分布

緑色の丸：結実確率に個体差があると仮定した場合の分布（平均値）

青色の領域： $\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の機能拡張やバグフィックスが続けられるとのことである（引用元）。

Kumada Seiya

仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

【C#】DataGridViewデータ表示速度比較

新卒の面接で逆質問を実施するときの恥ずかしさ

階層ベイズモデル

はじめに

ソースコード

問題設定

モデル化

個体差を考慮しない予測

個体差を考慮する予測

まとめ

Kumada Seiya

最近の記事

レイアウト解析

ChatGPTを用いた画像異常検知

クォータニオンと回転行列の関係

OpenCVで深層学習

OpenCVのcv::Matに関する小ネタ集

テンソルネットワークの入り口（その2）

ディープラーニングをやらないPyTorch入門

Bayes推定

時系列解析-1（MA過程）

次元の呪い

PandasをDaskで簡単に並列処理して高速化しよう

Learning Chainer with Examples 〜深層学習への入門〜

Google Vision APIでOCR

アーカイブ

カテゴリー