ディリクレ分布と多項分布
データ無しでStan、PyMC3を使う
ディリクレ分布は（１）式で与えらる。
多項分布は（２）式で与えられる。
ディリクレ分布には、確率変数
の和が１であるという制約があり、多項分布には確率変数
の和がNであるという同様の制約がある。このため、同時分布と周辺分布とはかなり様相が異なり得るが、このことをMCMCサンプリングによって調べてみた。ディリクレ分布の場合のStanスクリプト例とPyMC3スクリプト例、多項分布の場合のPyMC3スクリプト例を以下に示す。
ディリクレ分布（１）において、
のとき、
となり、
は一様分布である。しかし、周辺分布は、例えば
の値が小さいときは、他の確率変数の動ける領域は広くなり、の周辺確率は大きくなる。逆に、の値が大きいときは、他の確率変数の動ける領域は狭くなり、の周辺確率は小さくなる。すなわち、周辺分布は一様分布ではない。
式（１）においてK=5、
とおいた５変数の場合のStanスクリプトを用意して、周辺分布を調べるスクリプトをリスト１のように用意した。
リスト１．ディリクレ分布の周辺分布を調べるPyStan（Python用Stan）を用いたPythonスクリプト例。
import pystan
import matplotlib.pyplot as plt
import seaborn as sb
stan_code = """
parameters {
simplex[5] theta;
}
model {
theta ~ dirichlet([1.0, 1.0, 1.0, 1.0, 1.0]');
}
"""
sm = pystan.StanModel(model_code =
stan_code)
fit = sm.sampling()
plt.title('Marginal Distributions',
fontsize = 16)
for i in range(5):
sb.kdeplot(fit['theta'][:,i], clip = [0.0, 1.0], label = f'theta{i}')
plt.hist(fit['theta'][:,i],
density = True, alpha = 0.5)
plt.xlabel('theta', fontsize = 14)
plt.legend()
plt.savefig('StanDirichlet.png')
plt.show()
リスト１を実行すると、図１のグラフを得る。
図１ の周辺分布（Stan版）
同時分布では一様分布のが、周辺分布では左に偏った分布であることが分かる。
ディリクレ分布（１）において、のとき、
となり、は一様分布である。しかし、周辺分布は、例えばの値が小さいときは、他の確率変数の動ける領域は広くなり、の周辺確率は大きくなる。逆に、の値が大きいときは、他の確率変数の動ける領域は狭くなり、の周辺確率は小さくなる。すなわち、周辺分布は一様分布ではない。
式（１）においてK=5、とおいた５変数の場合のPyMC3スクリプトを用意して、周辺分布を調べるスクリプトをリスト２のように用意した。
リスト２ ディリクレ分布の周辺分布を調べるPyMC3を用いたPythonスクリプト例。
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sb
with pm.Model() as model_dirichlet:
theta =
pm.Dirichlet('theta', np.array([1.0, 1.0, 1.0, 1.0, 1.0]))
trace_theta =
pm.sample(3000)
plt.title('Marginal Distributions',
fontsize = 16)
for i in range(5):
sb.kdeplot(trace_theta['theta'][:,i], clip = [0.0, 1.0], label =
f'theta{i}')
plt.hist(trace_theta['theta'][:, i], density = True, alpha = 0.5)
plt.xlabel('theta', fontsize = 12)
plt.legend()
plt.savefig('DirichletPyMC3.png')
plt.show()
リスト２を実行すると、図２のグラフを得る。
図２ の周辺分布（PyMC3版）
同時分布では一様分布のが、周辺分布では左に偏った分布であることが分かる。
多項分布（式（２））において、すべての確率を等しく
と置くと、
となる。同時分布（３）は一様分布ではないが、
の和が一定Nであることから、確率変数の１つ、例えば
の値が小さいときは他の
は広い範囲の値を取ることができるが、の値が大きくなると他のの取り得る値の範囲は狭くなるというディリクレ分布と共通の性質がある。
カテゴリ数K=5、データ数N＝５、
と置いたときの確率変数の周辺分布を調べるスクリプトをリスト３のように用意した。
リスト３ 多項分布の周辺分布を調べるPyMC3を用いたスクリプト
import pymc3 as pm
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sb
with pm.Model() as check_pm:
x
= pm.Multinomial('x', n = 5, p = np.array([0.2, 0.2, 0.2, 0.2, 0.2]), shape =
5)
trace_x
= pm.sample(3000)
counts = np.zeros((5, 6))
for k in range(5):
x
= trace_x['x'][:,k]
for
i in range(len(x)):
counts[k][x[i]]
+= 1
plt.title('Marginal Frequencies',
fontsize = 14)
for k in range(5):
plt.plot(range(6),
counts[k], label = f'Category-{k}')
plt.xlabel('Coounts', fontsize = 12)
plt.ylabel('Frequency', fontsize = 12)
plt.legend()
plt.tight_layout()
plt.savefig('Multinomial.png')
plt.show()
リスト３を実行すると図３のグラフを得た。
図３ 多項分布の周辺分布
ディリクレ分布の場合と同様の傾向、すなわち、確率変数の値が小さいとその周辺確率は高いが、確率変数の値が大きくなると周辺確率は小さくなるという傾向が認められる。多項分布の場合は、式（３）に示されているように同時分布は一様分布ではないので、の値が０のときより１のときの方が周辺確率は高くなっていることが、サンプリングによるランダムな変動を除いて認められる。
なお、PyStaとPyMC3については＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞においてStan（第10、11、12章）およびPyMC3（第13、14、15章）で説明している。