平均値、スケールとｔ分布の自由度のベイズ推定
離散パラメータのMCMCによる事後分布推定
データの分布が正規分布に対して外れ値を含むと思われるとき、ｔ分布が用いられることがある。ｔ分布は、３つのパラメータ、位置パラメータmu、スケールパラメータsigma、自由度nu、をもつ。標準偏差は、nu>2のとき存在して
である。自由度は、正の実数であるが、ここでは１以上の整数値で考えて、離散パラメータとして扱う。
データが以下のように与えられたとする。
data = [185, 193, 209, 176, 211, 210, 208, 193, 195, 216,
178, 188, 209, 361, 190, 229, 246, 194, 162, 195,
215, 201, 217, 197, 218, 224, 209, 226, 189, 211,
213, 216, 215, 206, 233, 144, 186, 199, 225, 182,
209, 198, 191, 212, 177, 177, 146, 194, 188, 193]
このデータに対して、ｔ分布によるベイズ分析を行うPyMC3スクリプトを以下のように用意した。
with pm.Model() as t_model:
μ = pm.Normal('μ', mu = 0.0,
sigma = 1000.0)
σ = pm.HalfNormal('σ', sigma
= 1000.0)
v_df =
pm.DiscreteUniform('v_df', lower = 1, upper = 10)
t = pm.StudentT('t', nu =
v_df, mu = μ, sigma = σ, observed = data)
trace = pm.sample(3000)
分析を行う完全なスクリプトは、ウェブサイトの後ろにリスト１として用意した。なお、PyMC3の説明は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
リスト１を実行すると（Windowsの場合、Jupyter
Notebookで実行するとよい）、以下のようにMCMCサンプリングの結果が表示される。
図１.1
図1.2
図1.3
図1.4
上の分析結果とデータとの適合性を調べるために、データのヒストグラムと、パラメータを
と置いたときのｔ分布のグラフを描いてみた。プログラムは、リスト２に示すものである。
リスト２を実行すると、図２を得る。
図２
データと分析結果の適合性はよいと言えそうである。
import pymc3 as pm
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
data = [185, 193, 209, 176, 211, 210, 208, 193, 195, 216,
178, 188, 209, 361, 190, 229, 246, 194, 162, 195,
215, 201, 217, 197, 218, 224, 209, 226, 189, 211,
213, 216, 215, 206, 233, 144, 186, 199, 225, 182,
209, 198, 191, 212, 177, 177, 146, 194, 188,
193]
with pm.Model() as t_model:
μ = pm.Normal('μ', mu = 0.0,
sigma = 1000.0)
σ = pm.HalfNormal('σ', sigma
= 1000.0)
v_df =
pm.DiscreteUniform('v_df', lower = 1, upper = 10)
t = pm.StudentT('t', nu =
v_df, mu = μ, sigma = σ, observed = data)
trace = pm.sample(3000)
print(az.summary(trace))
# 日本語表記用の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
az.plot_posterior(trace['μ'])
plt.title('μの事後分布', fontsize = 20)
plt.savefig('Figure_mu.png')
plt.show()
az.plot_posterior(trace['σ'])
plt.title('σの事後分布', fontsize = 20)
plt.savefig('Figure_sgm.png')
plt.show()
plt.plot(range(len(trace['v_df'])),
trace['v_df'], 'b', alpha = 0.3)
plt.title('v_df のトレース', fontsize = 20)
plt.xlabel('試行数', fontsize = 16)
plt.ylabel('サンプル値', fontsize = 16)
plt.savefig('Figure_v_df_trace.png')
plt.show()
unique, counts =
np.unique(trace['v_df'], return_counts = True)
print(unique, counts)
Mod_df = unique[np.argmax(counts)]
print('Mod_df =', Mod_df)
plt.hist(trace['v_df'], bins
=np.arange(0.5, 10.6, 1))
plt.title('v_df の事後分布\nMode =
{}'.format(Mod_df), fontsize = 20)
plt.xlabel('v_df', fontsize = 16)
plt.yticks([])
plt.tight_layout()
plt.savefig('Figure_v_df_hist.png')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as ss
data = [185, 193, 209, 176, 211, 210, 208, 193, 195, 216,
178, 188, 209, 361, 190, 229, 246, 194, 162, 195,
215, 201, 217, 197, 218, 224, 209, 226, 189, 211,
213, 216, 215, 206, 233, 144, 186, 199, 225, 182,
209, 198, 191, 212, 177, 177, 146, 194, 188, 193]
plt.hist(data, bins = np.arange(100,
300, 20), density = True)
x = np.linspace(50, 350, 1000)
f_distr = ss.t(df = 3, loc = 201,
scale = 17)
y = f_distr.pdf(x)
plt.plot(x, y)
plt.savefig('Figure_fit.png')
plt.show()