Up

不定パラメータのベイズ推定

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

多項分布

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

と与えたとき、パラメータ数は2個である。

データは、N回の試行において、2カテゴリの頻度が回と回であるとき、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

 

が直線の近傍に集まっていることが判る。個々の値としては不定であるが、確率を与える重みの比率の関係が反映されている。

重みの値を横軸にとり、確率を縦軸の値としてサンプル値をプロットすると図3を得る。

図3

 

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

図4

 

重みの値に関わらず比0.3の近傍に位置していることが判る。

 

上記のモデルでは、不定なモデルの例として重みから確率を算出しているが、確率を直接ディリクレ分布で与えると、この不定性はなくなる。

また、上記モデルでは、尤度関数の形がそのまま事後確率分布に反映されるように事前分布として一様分布を設定している。何か事前情報がある場合は、それに対応した事前分布が設定される。また、確率をsoftmaxモデルで与える場合は、事前分布はsoftmaxモデルの各パラメータに対して設定されることになる。

 

StanPythonで使う説明は、岡本安晴「いまさら聞けないPythonでデータ分析―多変量解析,ベイズ統計分析(PyStan,PyMC)―」に用意している。

 

 

Up