MCMCサンプリングにおける不定性とスィッチング
ベイズ分析における対応の試み
モデルのパラメーターの不定性の問題は、最小2乗法による重回帰分析の場合はコントラストの推定可能性の問題として扱われている（岡本、2014、pp. 39-42）。ベイズ分析におけるモデルは重回帰モデル以外においても不定性の問題が起こりうるが、簡単な２項分布モデルの場合においてパラメータの不定性の問題を見てみる。取り上げた２項分布モデルの場合、パラメータ推定におけるＭＣＭＣサンプリングのスィッチングの問題も観察されたので、その対応も試みた。スクリプトファイルは、まとめて圧縮ファイルindefswfiles.zipとした。
次のモデルを設定する。
物体AとBがあり、それぞれ＋とーの２つの内部状態をとるものとする。
観測されるのは、AとBをまとめた状態ABで、
Aの状態が＋でBの状態が＋のとき、あるいはAの状態がーでBの状態がーのときに、ABの状態は＋と観測さる。
Aの状態が＋でBの状態がーのとき、あるいはAの事態がーでBの状態が＋のときは、ABの状態はーと観測される。
Aの状態とBの状態は独立であるとし、それぞれの確率は以下のように表されるものとする。
このとき、
となる。
今、ｎ回の観測においてABが＋であった回数がｋ回であったとする。２項分布が適用できるとすると、
と置くことができる。
パラメータpは、データに対してモデル（３）により求めることができるが、pを与えるパラメータpaとpbは式（２）からは一意には決まらない。独立変数としてpaとpbをとったときのpの値の等高線図を描くと図１のようになっている。
図１
図１のグラフは、リスト１のスクリプトで描かれたものである。
また、式（２）において、paおよびpbを(1-pa)および(1-pb)と交換しても式（２）の値は変わらない。このため、MCMCサンプリングにおいてスィッチングが起こる。スィッチング対策例も示した。
モデル（３）によるベイズ分析を行うためにリスト２およびリスト３に示すスクリプトを用意した。リスト２では、Stanスクリプトのビルドと実行を行っている。Pystan3は、Windows上では直接実行できないので、WSLによりインストールしたUbuntu上で実行した。WSLではグラフィック描画ができないので、Stanスクリプトの実行により生成したMCMCサンプルの描画は、リスト２の実行で生成したMCMCサンプルをファイルに保存して、Windows上で実行したリスト３のスクリプトにより保存したMCMCサンプルを読み込んで分析・描画した。リスト２のスクリプトファイルとリスト３のスクリプトファイルを同じフォルダ内に置いておけば、リスト２の実行に続けてリスト３の実行を行うことは容易である。
リスト２の実行後、リスト３を実行すると、図２のグラフが表示される。
図２
パラメータpaとpbの事後分布は、０から１まで一様に広がっているが、これら２つのパラメータから算出されるpの分布は0.7付近に集まっている。式（３）によりパラメータpは、データに対応するモデルのパラメータとして定まるが、式（２）で定まるパラメータpaとpbは不定であることに対応している。個々のパラメータが不定のとき、パラメータの組み合わせによっては推定可能であることは、重回帰分析では推定可能性の問題として論じられている（岡本、2014、pp.39-42）。この意味で、パラメータpは推定可能であるといえる。
パラメータpa、pbの事後分布の裾が、０より左側、および１より右側に広がっているのは、KDEにより密度関数が推定されているためである。図２のWindowを閉じると、図３のヒストグラムが表示される。
図３
ヒストグラムでは、確率の事後分布は０から１の範囲に収まっている。
図３のWindowを閉じると、図４の散布図が表示される。
図４
点（pa, pb）の分布が、右上の領域と左下の領域に分かれているのは、paおよびpbと(1-pa)および(1-pb)の間のスィッチングであると考えられる。スィッチングを避ける方法として、ベイズ分析の場合は事前分布に制約を置くことが考えられるが、例をこの後で示した。
MCMCサンプリングのトレースプロット
パラメータpaとpbのMCMCサンプリングのトレース図を描くために、リスト４のスクリプトを作成して、Ubuntu 22において実行した。
図５に示すトレースプロットが得られた。
図５
パラメータpaとpbのスィッチングは、図４を見ると、右上のグループと左下のグループの間で起こっていると思われるので、これを防ぐために事前分布を次のように設定して、paとpbの組を左下のグループに制限してみた。
pa ~
uniform(0, 0.5);
pb ~
uniform(0, 0.5);
リスト５のスクリプトを、Python3.12/miniconda/Ubuntu24/WSL/Windows11において実行後、リスト６のスクリプトをPython3.13/Windows11において実行すると、図６のグラフが表示された。
図６
パラメータpaとpbの事後分布が左側に制限されていることがわかる。
図６のWindowを閉じると、図７のグラフが表示される。
図７
パラメータpaとpbの事後分布は、事前分布で制限された範囲でほぼ一様に分布しているが、これは不定性の反映であると思われる。
図７のWindowを閉じると、図８の散布図が表示される。
図８
散布図からは、スィッチングが起こっているようには見えない。
リスト１ 等高線描画スクリプト（draw_contour.py）
import matplotlib.pyplot
as plt
import numpy
as np
import matplotlib.cm as cm
x = np.linspace(0.0,
1.0, 100)
y = np.linspace(0.0,
1.0, 100)
X, Y = np.meshgrid(x,
y)
Z = X*Y + (1-X)*(1-Y)
# p = pa * pb + (1 -
pa) * (1 - pb)
fig, ax = plt.subplots()
CS = ax.contour(X,
Y, Z)
plt.xlabel('pa', fontsize=14)
plt.ylabel('pb', fontsize=14)
ax.clabel(CS, fontsize=10)
ax.set_title('p=pa*pa+(1-pa)*(1-pb)', fontsize=16)
plt.show()
リスト２ MCMCサンプリングスクリプト（indef_switch_ba.py）
import stan
import arviz
as az
import pickle
import pandas as pd
StanCode = """
data {
int
k;
int
n;
}
parameters {
real<lower=0, upper=1> pa;
real<lower=0, upper=1> pb;
}
transformed parameters {
real<lower=0, upper=1> p;
p =
pa*pb + (1-pa)*(1-pb);
}
model {
pa ~
uniform(0, 1);
pb ~
uniform(0, 1);
k ~
binomial(n, p);
}
"""
sm = stan.build(StanCode, data = {'n':10, 'k':7})
fit = sm.sample()
infer_data = az.from_pystan(posterior=fit,
posterior_model=sm)
smry = az.summary(infer_data)
print(smry)
df = fit.to_frame()
print(df.describe().T)
with open('df.pkl',
'wb') as f:
pickle.dump(df, f)
リスト３ MCMCサンプルの分析（indef_switch_anal.py）
import pickle
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn as sb
with open('df.pkl',
'rb') as f:
df
= pickle.load(f)
sb.kdeplot(df["pa"],
c='r', label='pa')
sb.kdeplot(df['pb'],
c='g', label='pb')
sb.kdeplot(df['p'],
c='b', label='p')
plt.xlabel('')
plt.legend()
plt.show()
plt.figure(figsize=(5,6))
plt.subplot(211)
plt.hist(df['pa'])
plt.xlabel('pa', fontsize=14)
plt.subplot(212)
plt.hist(df['pb'])
plt.xlabel('pb', fontsize=14)
plt.tight_layout()
plt.show()
plt.figure(figsize=(7,7))
plt.plot(df['pa'], df['pb'], 'o', markerfacecolor='none')
plt.xlabel('pa', fontsize=14)
plt.ylabel('pb', fontsize=14)
plt.show()
リスト４ Linux用スクリプト（indef_switch_linux.py）
import stan
import arviz
as az
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn as sb
StanCode = """
data {
int
k;
int
n;
}
parameters {
real<lower=0, upper=1> pa;
real<lower=0, upper=1> pb;
}
transformed parameters {
real<lower=0, upper=1> p;
p =
pa*pb + (1-pa)*(1-pb);
}
model {
pa ~
uniform(0, 1);
pb ~
uniform(0, 1);
k ~
binomial(n, p);
}
"""
sm = stan.build(StanCode, data = {'n':10, 'k':7})
fit = sm.sample()
infer_data = az.from_pystan(posterior=fit,
posterior_model=sm)
smry = az.summary(infer_data)
print(smry)
az.plot_trace(infer_data,
var_names=['pa', 'pb'])
plt.show()
df = fit.to_frame()
print(df.describe().T)
sb.kdeplot(df["pa"],
c='r', label='pa')
sb.kdeplot(df['pb'],
c='g', label='pb')
sb.kdeplot(df['p'],
c='b', label='p')
plt.xlabel('')
plt.legend()
plt.show()
plt.figure(figsize=(5,6))
plt.subplot(211)
plt.hist(df['pa'])
plt.xlabel('pa', fontsize=14)
plt.subplot(212)
plt.hist(df['pb'])
plt.xlabel('pb', fontsize=14)
plt.tight_layout()
plt.show()
plt.figure(figsize=(7,7))
plt.plot(df['pa'], df['pb'], 'o', markerfacecolor='none')
plt.xlabel('pa', fontsize=14)
plt.ylabel('pb', fontsize=14)
plt.show()
リスト５ スィッチングを避ける事前分布の設定（indef_nswitch_ba.py）
import stan
import arviz
as az
import pickle
import pandas as pd
StanCode = """
data {
int
k;
int
n;
}
parameters {
real<lower=0, upper=1> pa;
real<lower=0, upper=1> pb;
}
transformed parameters {
real<lower=0, upper=1> p;
p =
pa*pb + (1-pa)*(1-pb);
}
model {
pa ~
uniform(0, 0.5);
pb ~
uniform(0, 0.5);
k ~
binomial(n, p);
}
"""
sm = stan.build(StanCode, data = {'n':10, 'k':7})
fit = sm.sample()
infer_data = az.from_pystan(posterior=fit,
posterior_model=sm)
smry = az.summary(infer_data)
print(smry)
df = fit.to_frame()
print(df.describe().T)
with open('df_nsw.pkl',
'wb') as f:
pickle.dump(df, f)
リスト６ リスト５のスクリプトと組み合わせたMCMCサンプル表示スクリプト（indef_nswitch_alan.py）
import pickle
import numpy
as np
import pandas as pd
import matplotlib.pyplot
as plt
import seaborn as sb
with open('df_nsw.pkl',
'rb') as f:
df
= pickle.load(f)
sb.kdeplot(df["pa"],
c='r', label='pa')
sb.kdeplot(df['pb'],
c='g', label='pb')
sb.kdeplot(df['p'],
c='b', label='p')
plt.xlabel('')
plt.legend()
plt.show()
plt.figure(figsize=(5,6))
plt.subplot(211)
plt.hist(df['pa'])
plt.xlabel('pa', fontsize=14)
plt.subplot(212)
plt.hist(df['pb'])
plt.xlabel('pb', fontsize=14)
plt.tight_layout()
plt.show()
plt.figure(figsize=(7,7))
plt.plot(df['pa'], df['pb'], 'o', markerfacecolor='none')
plt.xticks(np.arange(0,
1.001, 0.1))
plt.yticks(np.arange(0,
1.001, 0.1))
plt.xlabel('pa', fontsize=14)
plt.ylabel('pb', fontsize=14)
plt.show()