次のコードのように、重みの未知のベクトルが確率質量関数になるように制約されているため、ディリクレ分布としてモデル化されている pymc3 で線形回帰モデルを実装しています。
with pm.Model() as model:
#prior on precision of normal likelihood
tau = pm.Gamma('tau', alpha=1, beta=1)
phi = np.empty(ncountries, dtype=object)
y = np.empty((nyears-1, ncountries), dtype=object)
for icountry, country in enumerate(countries):
#prior Dirichlet allocation for each country
phi[icountry] = pm.Dirichlet('mix_{c}'.format(c=country),
np.roll(mix, icountry),
shape=ncountries)
for iyear, year in enumerate(years[1:]):
suffix = '_{y}-{c}'.format(y=year, c=country)
previous_pop = Xs[iyear, :]
#likelihood
y[iyear, icountry] = pm.Normal('obs' + suffix,
mu=pm.Deterministic(
'mu' + suffix,
dot(phi[icountry], previous_pop)),
tau=tau,
observed=Ys[iyear, icountry])
実行して事後をサンプリングした後:
start = pm.find_MAP()
step = pm.Metropolis()
nsteps = 1000
trace = pm.sample(nsteps, step, start=start)
ディリクレ変数のトレースを分析したところ、それらの値が 1 に加算されないことがわかりました (以下は例です)。
array([[ 0.01029745, 0.00627394, 0.00996922, ..., 1.83955829,
0.00962185, 0.01020659],
[ 0.01029745, 0.00627394, 0.00996922, ..., 1.83955829,
0.00962185, 0.01020659],
[ 0.01029745, 0.00627394, 0.00996922, ..., 1.83955829,
0.00962185, 0.01020659],
...,
[ 0.02050308, 0.01685555, 0.01976797, ..., 1.92278065,
0.03956622, 0.00473735],
[ 0.01993214, 0.01632033, 0.01994876, ..., 1.92487858,
0.04078728, 0.00481424],
[ 0.01900882, 0.01528191, 0.02100671, ..., 1.92485693,
0.0395159 , 0.00524575]])
私は theano 変数に精通しておらず、ディリクレ RV が pymc3 でどのように表現されているかを調べるのが難しいことがわかりました...何か間違ったことをしていますか、それともトレースで返された値を正規化して合計が 1 になるようにする必要がありますか?
クイック更新
関数pm.find_MAP()
は一種の勾配降下最適化を採用しているようです。これは、ディリクレ分布からのドローを表すベクトルが確率質量関数であるという事実から生じる制約を考慮していません (その値は正でなければならず、それらの合計は 1 でなければなりません)。この制約は、アルゴリズムのサンプリング段階でも明らかに適用されず、尤度分布の精度がゼロに向かってドリフトするため、収束の問題が発生します。