不定パラメータのベイズ推定
データの自由度よりパラメータの数が多いとき、パラメータ値は一意には定まらない。このモデルが同定(identifiable)できないとき、パラメータ値を求めるために極値探索法を用いても収束せず、例えば、最尤法では最尤推定値は定まらない。しかし、ベイズ分析の場合は、このモデルが同定できないことは事後分布に反映される形で表されるので、パラメータの同定(identifiability)条件が満たされなくても分析可能である。このことを2カテゴリの多項分布で確認したものを以下に示す。
多項分布

において、 の場合を考える。このとき、パラメータを
の場合を考える。このとき、パラメータを

と与えたとき、パラメータ数は2個である。
データは、N回の試行において、2カテゴリの頻度が 回と
回と 回であるとき、1個の制約
回であるとき、1個の制約

があるので、データの自由度は である。
である。
したがって、パラメータの個数2はデータの自由度1より大きく、パラメータ と
と の値は不定である。
の値は不定である。
上の確率モデルを次のStanスクリプトで表す。
data {
        int
f[2];
}
parameters {
        real
w1;
        real
w2;
}
transformed parameters {
        vector[2]
p;
        p[1]
= w1 / (w1 + w2);
        p[2]
= w2 / (w1 + w2);
}
model {
        w1 ~
uniform(0.0, 1.0);
        w2 ~
uniform(0.0, 1.0);
        f
~ multinomial(p);
}
上のStanスクリプトのファイル名をBinCat.stanとして利用するPythonスクリプトを以下のように用意した。
import numpy as np
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
Data = {'f': [30, 70]}
sm = pystan.StanModel(file =
'BinCat.stan')
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
fit.plot()
plt.show()
w1 = fit['w1']
w2 = fit['w2']
plt.plot(w1, w2, 'b.', alpha = 0.3)
plt.plot([0.0, 3/7], [0.0, 1.0], 'g-',
linewidth = 3,
        
label = 'w1:w2 = 3:7')
plt.xlabel('w1', fontsize = 20)
plt.ylabel('w2', fontsize = 20)
plt.legend(loc = 4, fontsize = 20)
plt.show()
w1w2 = np.array(w1) / (np.array(w1) +
np.array(w2))
plt.xlabel('w1', fontsize = 15)
plt.ylabel('w1/(w1+w2)', fontsize =
15)
plt.plot(w1, w1w2, 'b.', alpha = 0.3)
plt.show()
sb.kdeplot(w1, w1w2)
plt.xlabel('w1', fontsize = 15)
plt.ylabel('w1/(w1+w2)', fontsize =
15)
plt.show()
上のPythonスクリプトでは、 、
、 である。実行すると、以下の結果を得る。
である。実行すると、以下の結果を得る。
まず、サンプリングのトレースグラフと分布グラフは図1のようになる。パラメータの事後分布が得られていることが判る。

図1
点推定値が不定である2つのパラメータ と
と の同時分布は図2に示される。
の同時分布は図2に示される。

図2
点 が直線
が直線 の近傍に集まっていることが判る。個々の値としては不定であるが、確率を与える重みの比率の関係が反映されている。
の近傍に集まっていることが判る。個々の値としては不定であるが、確率を与える重みの比率の関係が反映されている。
重み の値を横軸にとり、確率
の値を横軸にとり、確率 を縦軸の値としてサンプル値をプロットすると図3を得る。
を縦軸の値としてサンプル値をプロットすると図3を得る。

図3
図3の散布図を、カーネル密度関数で滑らかな曲面で表し等高線を描くと図4を得る。

図4
重み の値に関わらず比
の値に関わらず比 は0.3の近傍に位置していることが判る。
は0.3の近傍に位置していることが判る。
上記のモデルでは、不定なモデルの例として重みから確率を算出しているが、確率を直接ディリクレ分布で与えると、この不定性はなくなる。
また、上記モデルでは、尤度関数の形がそのまま事後確率分布に反映されるように事前分布として一様分布を設定している。何か事前情報がある場合は、それに対応した事前分布が設定される。また、確率をsoftmaxモデルで与える場合は、事前分布はsoftmaxモデルの各パラメータに対して設定されることになる。
StanをPythonで使う説明は、岡本安晴「いまさら聞けないPythonでデータ分析―多変量解析,ベイズ統計分析(PyStan,PyMC)―」に用意している。