書記が数学やるだけ#781 正規分布モデルのベイズ統計

鈴華書記（Writer Rinka）

2023年10月20日 08:42

続いて正規分布モデル，パラメータが2個に増えた場合について。

問題

説明

分散既知の場合は平均のみを推定する問題に帰着され，平均の共役事前分布として正規分布を設定する。

実際には分散未知の場合を考えることが多い，ここで分散の共役事前分布は逆ガンマ分布である。

逆ガンマ分布の確率密度関数を以下に示す。

解答

分散既知のモデルにおいて事後分布を求める。ここで分散の逆数を精度といい，今後もよく用いられる。

次に予測分布を求める。

分散既知の場合について実際に計算してみる。

分散未知の場合について複雑な計算は省略するが，確かに事後分布が正規分布となっている。

分散未知の場合について計算してみる。

では実装例を見ていく。今回はテストの点数を正規分布で近似する問題である。

import pymc3 as pm
import arviz as az
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn-darkgrid')

大体は50点付近で，1人だけ100点という外れ値がいる場合を考える。

#標本
data = np.array([50,45,55,60,40,100])
#カーネル密度推定
sns.kdeplot(data)

ではモデルを構築していく。ここで分散の事前分布に半正規分布を用いているのは，Osvaldo Martin『Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド』に合わせている。

またAndrew Gelman(2006)"Prior distributions for variance parameters in hierarchical models"では，逆ガンマ分布でなく半コーシー分布を使うのが良いことが示されている。さらに，"Prior Choice Recommendations"では事前分布の選び方がまとまっている。

#正規分布モデル
with pm.Model() as model_g:
#平均：一様分布
mu = pm.Uniform('mu', 40, 100)
#標準偏差：半正規分布
sigma = pm.HalfNormal('sigma', sd=10)
#標本
y = pm.Normal('y', mu=mu, sd=sigma, observed=data)

trace_g = pm.sample(1100)
chain_g = trace_g[100:] pm.traceplot(chain_g)
plt.figure()

結果をまとめておく。

az.summary(chain_g)

事後分布は以下の通りで，平均が外れ値の分だけ引きずられているのがわかる。

az.plot_posterior(chain_g)

続いて正規分布の代わりにt分布を採用する，これにより頑健推定を行うことができる。t分布では平均・分散のほかに自由度を指定する。

この自由度が小さい場合は裾の広い分布となり，自由度が十分大きい場合には正規分布に近似できる。

#t分布モデル
with pm.Model() as model_t:
#平均：一様分布
mu = pm.Uniform('mu', 40, 100)
#標準偏差：半正規分布
sigma = pm.HalfNormal('sigma', sd=10)
#自由度：指数分布
nu = pm.Exponential('nu', 1)
#標本
y = pm.StudentT('y', mu=mu, sd=sigma, nu=nu, observed=data)

trace_t = pm.sample(1100)
chain_t = trace_t[100:]
pm.traceplot(chain_t)
plt.figure()