十の並列した脳

何でも勉強する,毎週月木曜に投稿予定

「Pythonによるベイズ統計モデリング」まとめ #4

Pythonによるベイズ統計モデリングPyMCでのデータ分析実践ガイド― 」:

 

Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド

Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド

 

前回↓ 

 

ryosuke-okubo.hatenablog.com

 

 

#4では第6章と第7章を扱う。

 

 

6章 モデル比較

過剰適合と過少適合

次の図は,同じデータに対して線形モデルと非線形モデルを当てはめたものである。

f:id:ryosuke_okubo:20190212144411p:plain

果たしてどちらが,よりよいモデルだろうか?一見すると,右の方がデータへの当てはまりがよい。しかし,これは背後にある分布をとらえられているのか?新たにデータを加えてみると,様相が変わる。

f:id:ryosuke_okubo:20190212145048p:plain

まあ極端な例だけれども。

青点が新たに加えられたデータであり,右のモデルはそれについて当てはまりがよいとは言えない。一方で左のモデルは単純ではあるが,新たに加えられたデータに対しても当てはまりがよい。このように,モデルが複雑すぎると過剰適合を起こし,未知のデータに対する予測能が悪くなる。しばしば機械学習の文脈で過学習と言われるものである。

かといって単純すぎるのもよくない。上の例で,0次式でモデルをつくると,かえってデータの特徴がつかめなくなる。これを過少適合という。

 

まとめると,モデルの複雑さ予測能とではトレードオフの関係にあるといえる(補足:その複雑さに意味があるのかを考えるのも大事)。

f:id:ryosuke_okubo:20190212150717p:plain

 

事前分布の正則化

過剰適合を回避する1つの方法として,事前分布を正則化する。

例:

  • チコノフ正則化
  • リッジ回帰
  • ラッソ回帰

 

予測精度

モデル評価には以下の方法がある。

  • 交差確認:データを分割して,片方をモデル作成に,もう片方をモデル評価に用いる。
  • 情報量基準:モデル評価の近似式。
    • 逸脱度
    • 赤池情報量基準(AIC
    • 逸脱度情報量基準(DIC)
    • 広く使える情報量基準(WAIC)
    • パレート平滑化重点サンプリング(PSIS)と1個抜き交差確認(LOOCV)
    • ベイジアン情報量基準(BIC

 

ベイズファクター(BF)

{\displaystyle BF = \frac{p(y|M_0)}{p(y|M_1)}} ,ここでM_0,M_1はそれぞれ比較対象とするモデルを示す。ベイズファクターの意味として,BF \gt 1のときは,M_0M_1よりもデータをうまく説明していることを示している。情報量基準と似た概念と思ってもよい。

 

ベイズファクターを計算するためには,

{\displaystyle \frac{p(y|M_0)}{p(y|M_1)} =\frac{p(M_0|y)p(M_1)}{p(M_1|y)p(M_0)}}

とする(ベイズの定理より)。

 

コイン投げを例に説明する。100回投げて30回表が出る事例において,

  • モデル0:Beta(4,8),図の青線
  • モデル1:Beta(8,4),図のオレンジ線

のどちらがモデルとして優れているだろうか?

f:id:ryosuke_okubo:20190212205143p:plain

 コード例

mport pymc3 as pm
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')

#100回投げて30回表が出る
coins = 100
heads = 30
y_d = np.repeat([0, 1], [coins-heads, heads])

with pm.Model() as model_BF:
    p = np.array([0.5, 0.5])
    model_index = pm.Categorical('model_index', p=p)
    m_0 = (4, 8) #モデル0,Beta(4,8)
    m_1 = (8, 4) #モデル1,Beta(8,4)
    m = pm.math.switch(pm.math.eq(model_index, 0), m_0, m_1)
  
    theta = pm.Beta('theta', m[0], m[1])
    y = pm.Bernoulli('y', theta, observed=y_d)
    trace_BF = pm.sample(5500, njobs=1)

chain_BF = trace_BF[500:]
pm.traceplot(chain_BF)

plt.figure()

f:id:ryosuke_okubo:20190212205517p:plain

\thetaの事後分布をみるに,モデル0の方が近そうである。

ではベイズファクターを求めて,数量的なモデル評価をしてみる。

pM1 = chain_BF['model_index'].mean()
pM0 = 1 - pM1
BF = (pM0/pM1)*(p[1]/p[0])
print("pM0 =",pM0) #pM0 = 0.954
print("pM1 =",pM1) #pM1 = 0.046
print("BF =",BF) #BF = 20.73913043478261

 ベイズファクターの値は約21となった。これでモデル0がモデル1より優れていることが数量的に示すことができた。


7章 混合モデル 

これまでの章の応用例なので,概要のみ説明する。

 

今まで確率分布として正規分布や2項分布を用いてきたが,実際の現象の中にはより複雑な分布をとるものがある。そのような現象は,単純な分布を組み合わせた混合モデルを仮定することで説明できるものがある。

例えば成人の身長の分布は二峰性の分布をとる。これは男性グループと女性グループの組み合わせにより説明できる。

混合モデルを記述するためのモデルとして

などが用いられる。

 

次回↓

作成中