ベイズ統計学の概論的紹介

ベイズ統計学の 概論的紹介 東京工業大学 数理･計算科学系 林 直輝 本スライド や資料の二次配布は常識の範囲でお願いし ます。 1

本資料の位置づけ • 本資料は2019/08/06に東京工業大学における実習「数理・計算科学 チュートリアル実践第一」[1]で使用した講義資料を追記修正したもの である。 • 上記実習は修士課程学生向けのチュートリアル講義を行うものであり、 本資料とほぼ同一なものが[2]において公開済みである。 [1]http://www.ocw.titech.ac.jp/index.php?module=General&action=T0300&GakubuCD=4&GakkaCD=34222 2&KeiCD=22&course=22&KougiCD=201921140&Nendo=2019&lang=JA&vid=03 [2]https://www.slideshare.net/naokihayashi714/ss-161469671 次ページ以降が講義資料である 2

• http://www.ocw.titech.ac.jp/index.php?module=General&action=T0300&GakubuCD=4&GakkaCD=342222&KeiCD=22&course=22&KougiCD=201921140&Nendo=2019&lang=JA&vid=03
• https://www.slideshare.net/naokihayashi714/ss-161469671

本チュートリアルの事前知識と概要 • 教養レベルの解析学と線型代数学と確率論と統計学を仮定 • ベイズ推論の基礎概念から理論研究の成果までを概論的に紹介 ‒ ベイズ推論のフレームワーク（パラメータの分布を使って予測分布を導く） ‒ 代表的な推論アルゴリズム（ﾒﾄﾛﾎﾟﾘｽ=ﾍｲｽﾃｨﾝｸﾞ法、ｷﾞﾌﾞｽｻﾝﾌﾟﾘﾝｸﾞ） ‒ ベイズ統計理論（汎化損失と自由エネルギーの理論とモデルの評価方法） • 資料中で参照するJupyterNotebookは下記参照 http://nhayashi.main.jp/codes/BayesStatAbstIntro.zip https://github.com/chijan-nh/BayesStatAbstIntro 概論的：証明は扱わず結果のみを紹介…… 3

• http://nhayashi.main.jp/codes/BayesStatAbstIntro.zip
• https://github.com/chijan-nh/BayesStatAbstIntro

ベイズ推論の フレームワーク 4

ベイズ推論のフレームワーク • ベイズ推論の流れは次のようになる： ‒ ステップ1 モデルと事前分布を定義 ベイズ推論とは 「データを発生している真の分布は予測分布だろう」 と推論すること ◼ モデル：データの確率分布。人間が作る。パラメータを持つ ◼ 事前分布：モデルが持つパラメータの確率分布。人間が作る ‒ ステップ2 事後分布を計算 ◼ 事後分布：データを観測した後のパラメータの分布。アルゴリズムが作る ‒ ステップ3 予測分布を計算 ◼ 予測分布：データを観測した後の未知データも含むデータの分布。アルゴリズムが作る 5

事前分布とは何か？ • 事前分布とは何か？ 事前知識？ 主観？ モデルを使ってデータ の分布を推論します • 人間が作ったモデルのパラメータの分布？？？？ モデル データ 人間 6

事前分布とは • 事前分布とは、データをモデルが観測する前のパラメータの分布 ‒ データだけでなくパラメータも確率変数であるというモデリングをする ‒ 事前分布も統計モデルの一部であるというだけ データ モデル𝑝 𝑥 𝑤 モデルを使ってデータの 分布を推論します。 事前分布は𝜑 𝑤 を使っ てみます。 人間 7

尤度とは • 最尤推定で最大化していた目的函数（尤度函数） • ベイズ推論でも使う • 単に、各データ点でのモデルの値の総積である どんなにデータ数が多くても、 尤度が正規分布で近似でき るとは限りません。 最大化点が無数あったり存在 しないこともあります。 8

事後分布とは • 事後分布とは、データをモデルが観測した後のパラメータの分布 ‒ モデルがデータを観測する確率→尤度 ‒ 観測前のパラメータの確率 →事前分布 モデル𝑝 𝑥 𝑤 事前分布𝜑 𝑤 データ 尤度𝑃 𝑋 𝑛 𝑤 事後分布 𝑃 𝑋 𝑛 |𝑤 𝜑 𝑤 𝑛 𝜓 𝑤 𝑋 := 𝑍 𝑋𝑛 9

予測分布とは • 予測分布とは、データをモデルが観測した後のデータの分布 ‒ データを観測したことでパラメータの分布が事後分布に更新された ‒ この事後分布を使ってモデルを更新する データ モデル𝑝 𝑥 𝑤 事前分布𝜑 𝑤 尤度𝑃 𝑋 𝑛 𝑤 事後分布 𝑛 𝑃 𝑋 |𝑤 𝜑 𝑤 𝑛 𝜓 𝑤 𝑋 := 𝑍 𝑋𝑛 予測分布 𝑝 ∗ 𝑥 𝑋 𝑛 ≔ න𝑝 𝑥 𝑤 𝜓 𝑤 𝑋 𝑛 𝑑𝑤 モデル𝑝 𝑥 𝑤 事後分布𝜓 𝑤 𝑋 𝑛 モデルを事後分布で平均したもの → データの情報を捨てずに予測する ※最尤推定ではデータの情報を 最尤推定量という1点に潰してしまう 10

ベイズ推論のフレームワーク（再掲） • ベイズ推論の流れは次のようになる： ‒ ステップ1 モデルと事前分布を定義 ◼ モデル：データの確率分布。人間が作る。パラメータを持つ ◼ 事前分布：モデルが持つパラメータの確率分布。人間が作る ‒ ステップ2 事後分布を計算 ◼ 事後分布：データを観測した後のパラメータの分布。アルゴリズムが作る ‒ ステップ3 予測分布を計算 ◼ 予測分布：データを観測した後の未知データも含むデータの分布。アルゴリズムが作る 11

ベイズ推論の代表的なタスク • メーカー・工場製品の品質管理 • ロケットやロボットの制御・自己位置推定 • 新薬や新素材開発のための物質構造解析 • 推薦システムの構築 ※これは筆者とは別人です • セマンティック検索の実現 • 広告効果モデル・マーケットリサーチの実施 • （医学・生物学・地学・心理学といった分野でも盛ん） 12

代表的なアルゴリズム 13

事後分布の実現 ベイズ更新 データが次々得られて、何度 も更新したいときは特に！ 懸念 ベイズ更新により分布の形が複雑になってしまうと、ベイズ 更新が行えなくなったり、行いにくくなったりしそう。 14

事後分布の実現 ෠ 𝜑 𝑤 = 𝑓 𝑤 𝜙 , 𝜓 𝑤 = 𝑓 𝑤 𝜙)と書ける密度函数𝑓が存在 共役事前分布 ベイズ更新𝜑 ↦ 𝜓にあたり、分布の種類が変わらないような 事前分布を共役事前分布と呼ぶ。 15

共役事前分布の例 このとき，予測分布も解析的に計算可能 【引用元】全人類がわかる統計学 https://to-kei.net/bayes/conjugate-prior-distribution/ 16

• https://to-kei.net/bayes/conjugate-prior-distribution/

より一般に…… • 実問題では共役事前分布が存在するとは限らないモデルを使う • 今日広く使われる「構造を持つモデル」に対しては、いずれも直接 の共役事前分布はない ‒ 混合分布、トピックモデル、階層ベイズ…… • 単純なケースでも、モデルを正規分布、事前分布を一様分布、 とすることもあるが、共役事前分布にならない • そもそも事前分布は自由に設計できるはず • 自由に設計するためには、共役とは限らない場合においても利用 可能な計算手法（と数理科学）が必要である 17

一般の事後分布の実現 • 共役じゃないとしても事後分布は定義可能 • 定義通りに計算？ → 難しい この積分はいつでも 解けるとは限らない 18

事後分布の近似 • 一般に事後分布を厳密に計算することは難しい • 近似計算は？ 目の前にはコンピュータが…… 事後分布の近似方法の一つ、 マルコフ連鎖モンテカルロ法 19

MCMCの概要 • マルコフ連鎖モンテカルロ法（Markov Chain Monte-Carlo method, MCMC）とは？ • 求める確率分布を均衡分布として持つマルコフ連鎖を作成することをも とに、確率分布のサンプリングを行うアルゴリズムの総称[Wikipedia, 2019/10/30閲覧] ？？？？ 20

MCMCの概要 • 用語を整理 ‒ マルコフ連鎖：前の状態にしか依存しない時系列。 簡単な天気予報や文書生成など ‒ モンテカルロ法：いわゆるモンテカルロシミュレーション。ある範囲の乱数を大 量に生成することでその範囲の面積や体積を計算する。円周率πのナイーヴ な計算など 半径1の円の面積はπ せやな わかる 知らんけど 次の発話は今の発話に依存するが、 それより前の発話には依存しないと仮定 一様乱数を生成して円内に 入った個数の割合が面積比に 近づいていく →πに近づいていく 21

MCMCの概要 • MCMCはマルコフ連鎖をモンテカルロ法で作る • 何がうれしいの？ • 実はマルコフ連鎖はある確率分布（定常分布）に収束する → 事後分布に収束するようなマルコフ連鎖を作ればよい • 「事後分布に収束」させるマルコフ連鎖を作る方法がMCMCである！ ‒ 無限に続ければ収束する保証がある 22

MCMCの概要 • MCMCを使うと確率分布からサンプリングできることはわかった • 目的の確率分布からサンプリングするMCMCのアルゴリズムは何 だろう？ • 本講座では次の2つを紹介 ‒ メトロポリスヘイスティング法 ‒ ギブスサンプリング 23

メトロポリスヘイスティング法 • メトロポリスヘイスティング (Metropolis Hasting, MH) 法 • 汎用的なMCMCアルゴリズムの1つ，最も基本的なもの 24

メトロポリスヘイスティング法の原理 • MCMCの目的：ある確率分布に従う確率変数の列の作成 • ある分布に従うということは，局所的に従っていてもダメ 目的の分布から満遍なくサンプリングする必要がある！ 25

メトロポリスヘイスティング法の原理 • MCMCの目的：ある確率分布に従う確率変数の列の作成 • ある分布に従うということは，局所的に従っていてもダメ 目的の分布から満遍なくサンプリングする必要がある！ • 目的の分布が一様とは限らないので，適切に設計しなければならない • そのために，事後分布の「ポテンシャルエネルギー」を考慮する 26

メトロポリスヘイスティング法の原理 • 満遍なくサンプリングするには？ → 程よいステップ幅εで分布の「エネルギー面」を動き続ける 1. 2. 3. 4. 初期位置wを決める w’ = w + d, d～N(0,ε^2) として次の位置の候補w’を求める 「ある基準」で次の位置をw’にし，そうでなければwに留まる 2に戻る ここでエネルギーが役に立つ 27

メトロポリスヘイスティング法の原理 • 「ある基準」とは ‒ 疑問：動く必要があるなら常に動き続ければいいのでは？ ‒ 全くのランダムウォークではMCMCにならない（一様分布になる） ‒ 分布の形状（エネルギー面）を考慮しないと 事後分布からのサンプリングにならない ‒ 現在位置と次の位置候補のエネルギー差H(w’)-H(w)を計算する ◼ エネルギー：𝜓 𝑤 𝑋 𝑛 ∝ exp −𝐻 𝑤 なる函数H(w) ‒ 差が負のときは必ず移動する ‒ 差が非負のときは確率exp(-差)で移動する ◼ 差が非負のときにも動かないと， 停留点から動かなくなってしまい サンプリング箇所が偏ってしまう 28

メトロポリスヘイスティング法のアルゴリズム • 準備： ‒ エネルギー：𝜓 𝑤 𝑋 𝑛 ∝ exp −𝐻 𝑤 なる函数H(w)を用意する ‒ ステップ幅εを固定する • サンプリング： 1. パラメータの初期値wを決める 2. w’ = w + d, d ～ N(0,ε^2) 3. 確率min 1, exp −𝐻 𝑤 ′ + 𝐻 𝑤 そうでなければ留まる 4. 2に戻る 5. 以上を所定の回数繰り返す でw←w’（移動）とし， 例として，MCMC_Converg.ipynb 29

メトロポリスヘイスティング法の注意点 • ステップ幅を事前に決めておく必要がある • 実際にはいくつかの値を検討してサンプリングを成功させる必要がある • 成功させる規準として，繰り返し回数のうち移動した回数が2割から3割を 占めるように調整させるヒューリスティクスがある（？？） ‒ この割合を採択確率という ‒ 理論的には採択確率はいくつでもよい（無限列で収束する） • より高度なMHとして，No U Turn Sampling (NUTS) 法があり，エネルギーの 微分可能性を仮定するが、こちらはステップ幅を自動的に決定可能 ‒ 確率的プログラミング言語Stanを使うと簡単にNUTSを使ったMCMCが可能 • まとめ：適用可能な状況は多いが，細かい設定が必要 ‒ エネルギー函数Hは尤度と事前分布の積に-logをつければすぐ求まる 30

メトロポリスヘイスティング法の注意点 • 無事ステップ幅を決定しても， 実は先ほどのサンプリングには問題が残っている • １つは初期値の影響である． マルコフ連鎖の定常分布は本来初期値に影響されないが， 有限のサンプリングを行う実際の状況では影響を受けてしまう ‒ このためにバーンイン (burn-in) という処理がされる． サンプリング初期からいくつかのサンプルは用いないというものである． ‒ MCMC_Converg.ipynbでは最初の2万サンプルを捨てている 31

メトロポリスヘイスティング法の注意点 • 無事ステップ幅を決定しても， 実は先ほどのサンプリングには問題が残っている • ２つ目は相関性である．サンプリング時， 新たなサンプルが前のサンプルの位置などの影響を受けており， 強すぎる系列相関になってしまう． ‒ このためにシニング (thining) という処理がされる． サンプリングを一定の周期で行うというものである． ‒ MCMC_Converg.ipynbでは周期を20としている 32

ギブスサンプリング • MH法にはステップ幅の問題があった • また，エネルギー函数を追跡してサンプリングする必要があった • 他にサンプリング手法はないだろうか？ → ギブスサンプリング 33

ギブスサンプリング • ギブスサンプリング (Gibbs Sampling, GS) もMCMCの1つ • MH法でエネルギーの追跡などイロイロな処理が出てきたが， これは事後分布が部分的にも簡単に書くことができないためである • 事後分布そのものはわからなくとも， 各変数については解析的な分布がわかる場合はGSが有効である この「簡単になった」分布からのサンプリングを基にして事後分布 に従うサンプルを得るのがGS。解析的に計算しやすくするために 共役事前分布を意識したモデリングをすることが多い 34

ギブスサンプリング • w1及びw2という2つのパラメータがある • 事後分布𝜓(𝑤1 , 𝑤2 |𝑋 𝑛 ) は解析的に分からなくとも， それぞれの条件付き確率分布𝜓1 (𝑤1 |𝑤2 , 𝑋 𝑛 )と𝜓2 (𝑤2 |𝑤1 , 𝑋 𝑛 )は 簡単な分布で書けるとする ‒ 例えば正規分布であるなど（右図） • アルゴリズムは極めて簡明である： 1. 𝑤1 の初期値を決める 2. 𝑤2 ~ 𝜓2 (𝑤2|𝑤1, 𝑋 𝑛 ) 3. 𝑤1 ~ 𝜓1(𝑤1 |𝑤2 , 𝑋 𝑛 ) w1 w2 2つの分布から乱数生成 を繰り返すだけ！！ 4. 2と3を所定回数繰り返す 35

ギブスサンプリングの注意点 • ギブスサンプリングは導出できれば非常に簡明 • 導出できるのはある程度限られた状況になる ‒ 直接の共役事前分布ではないが， その性質が利用できるように帰着できる事前分布を使うことになる ‒ 混合正規分布やトピックモデルといった実際的な構造を持つモデルでは， 隠れ変数を導入してGSの導出が可能になることがある • 尤度と事前分布からすぐ計算できるMHのエネルギー函数とは異なり， GSの条件付き分布の導出は数式的に非自明 • MH同様にバーンインとシニングは必要 • まとめ：使える状況は限られるが使えるときは非常に強力 ‒ 簡単な分布からの乱数生成を指定回数繰り返すのみである 36

MCMCと予測分布 • MCMCより，事後分布𝜓 𝑤 𝑋 𝑛 からサンプル 𝑤𝑘 𝐾 𝑘=1 が得られる • 𝐾 → ∞のとき，函数𝑓 𝑤 についてほとんど至る所 • MCMCを用いた時の予測分布𝑝∗ 𝑥 は数値的には下記を用いる： 1 𝐾 ‒ 𝑝∗ 𝑥 と σ𝐾 𝑘=1 𝑝 𝑥 𝑤𝑘 が𝐾や 𝑤𝑘 によってどれほど異なるかは一般には未解明： 37

実験の前に…… • 今回取り上げるモデルの説明 38

混合ポアソン分布 • たいていの有名かつ簡単な確率分布は1つの峰を持つか平坦 • 複数峰があるデータにはどう対応するのだろうか？ 39

混合ポアソン分布 • たいていの有名かつ簡単な確率分布は1つの峰を持つか平坦 • 複数峰があるデータにはどう対応するのだろうか？ 複数の分布を 重ね合わせる 40

混合ポアソン分布 • 数式モデルは ‒ 𝜆 = 𝜆1 , … , 𝜆𝐾 :1つ1つのポアソン分布の平均 ‒ 𝜋 = 𝜋1 , … , 𝜋𝐾 :混合比、すべて𝜋𝐾 ≥ 0かつσ𝐾 𝑘=1 𝜋𝑘 = 1 • この2つのパラメータの事前分布を考える • メトロポリスヘイスティングでは自由に • ギブスサンプリングでは ‒ 𝜆: ガンマ分布 ‒ 𝜋: ディリクレ分布 混合ポアソン分布は構造を持つモデル （特異モデル）ですので、最尤法や事後 確率最大化法に比べてベイズ推論は 精度が良い手法になります 41

混合ポアソン分布のギブスサンプリング • 混合ポアソン分布のGSを導出するには、更に隠れ変数を定義します • 隠れ変数s_ik = i番目のサンプルがk番目のコンポーネントに属すとき 1、そうでなければ0 隠れ変数 混合比 データ 各コンポーネントの 平均パラメータ 42

混合ポアソン分布のギブスサンプリング • 隠れ変数を導入すると、数式としては ‒ 隠れ変数𝑠はk番目が1で他が0のK次元ベクトル（K次元onehotベクトルという） ‒ すべての𝐾次元onehotベクトル𝑠について周辺化すると、元の𝑝 𝑥 𝜆, 𝜋 に戻る • ギブスサンプリングを導出するにはこの隠れ変数の導入が必要 43

混合ポアソン分布のギブスサンプリング • 導出の概略として、 ‒ 隠れ変数が与えられたときのパラメータの分布𝜓𝜆 𝜆 𝑆 𝑛 , 𝑋 𝑛 , 𝜓𝜋 𝜋 𝑆 𝑛 , 𝑋 𝑛 ‒ パラメータが与えられたときの隠れ変数の分布𝜓𝑆 𝑆 𝑛 |𝜆, 𝜋, 𝑋 𝑛 ‒ これら2つの条件付き確率分布が共役事前分布を用いて簡単に書けることを使う 隠れ変数 混合比 データ 各コンポーネントの 平均パラメータ 44

混合ポアソン分布のギブスサンプリング （導出は省略するが）結果として、次のアルゴリズムが得られる 1. パラメータ 𝜆, 𝜋 に初期値を設定 2. 𝜓𝑆 𝑆 𝑛 |𝜆, 𝜋, 𝑋 𝑛 は多項分布であり、そこからサンプリング 3. 𝜓𝜆 𝜆 𝑆 𝑛 , 𝑋 𝑛 はガンマ分布であり、そこからサンプリング 4. 𝜓𝜋 𝜋 𝑆 𝑛 , 𝑋 𝑛 はディリクレ分布であり、そこからサンプリング 5. サンプリングした 𝜆, 𝜋 を用いて2.に戻る 6. 以上を所定回数繰り返す 45

実験 • 構造を持つモデルでMHとGSをそれぞれ動かしてみる • モデル：混合ポアソン分布 (Poisson Mixture Model, PMM) • MH： PMM_MH.ipynb • GS： PMM_GS.ipynb 46

ベイズ統計理論 47

推論して得られた結果 • ここまで、様々なベイズ推論の方法を扱ってきた • 推論すれば何か結果は得られる • ……so what? • 得られた推論結果を導き出したモデルの評価を行いたい 48

推論して得られた結果 • ベイズ推論に限らず、統計的推論は正しくない。 • 推論の間違いは「正しく」知ることができる！ ‒ どうやって？ 真の分布（実世界）はわからないのに…… → 数理科学の出番 ‒ 真の分布が分からなくても、データとモデルから推論の間違いを知ることができる ◼ 真の分布がどんなものであっても成り立つ普遍法則を作ればよい ‒ 例えば ◼ ベイズ推論では、真の分布が分からなくてもデータとモデルから計算できる量を使って、推論 結果の間違いがどのような挙動を持つか導出することができる ◦ 最尤推定でも似たような定理が証明できるが、モデルと真の分布に強い制約がつく 49

推論して得られた結果 • ベイズ推論に限らず、統計的推論は正しくない。 • 推論の間違いは「正しく」知ることができる！ ‒ どうやって？ 真の分布（実世界）はわからないのに…… これにまつわる → 数理科学の出番 数理科学の成果を紹介 ‒ 真の分布が分からなくても、データとモデルから推論の間違いを知ることができる 定理が成立する仮定の下で正しい ◼ 真の分布がどんなものであっても成り立つ普遍法則を作ればよい ‒ 例えば ◼ ベイズ推論では、真の分布が分からなくてもデータとモデルから計算できる量を使って推論結 果の間違いがどのような挙動を持つか計算する理論を導出することができる ◦ 最尤推定でも似たような定理が証明できるが、モデルと真の分布に強い制約がつく 50

記号の復習 • モデル：𝑝 𝑥 𝑤 • 事前分布：𝜑 𝑤 • データ：𝑋 𝑛 = 𝑋1 , … , 𝑋𝑛 • 尤度：𝑃 𝑋 𝑛 |𝑤 • 事後分布：𝜓 𝑤 𝑋 𝑛 • 周辺尤度、エビデンス：𝑍 𝑋 𝑛 • 予測分布：𝑝 ∗ 𝑥 given data を強調して𝑝∗ 𝑥 𝑋 𝑛 とも書く 51

準備：KL情報量 • 最も一般的な分布間の非類似度としてカルバック・ライブラ (KL) 情報量がある： 𝑞 𝑥 𝐾𝐿 𝑞||𝑝 = න 𝑞 𝑥 log 𝑑𝑥 𝑝 𝑥 ‒ という定義式で、qからpへのKL情報量という • 情報理論や統計力学の負の相対エントロピーと同値 ‒ 情報理論では情報源と受信機の非類似度 ‒ 統計学・機械学習ではデータの分布とモデルの非類似度 • どちらも送信と受信で非対称な関係なため、KLは非可換𝐾𝐿 𝑞||𝑝 ≠ 𝐾𝐿(𝑝| 𝑞 ‒ 距離にはならないが なる正定値性を持つ。 𝐾𝐿 𝑞||𝑝 ≥ 0, 𝐾𝐿(𝑞| 𝑝 = 0 ⇔ 𝑝 = 𝑞 52

統計的推論の目的と定式化 • 「正確な予測」（予測重視） ‒ データに基づいて未知データの振る舞いを（可能な限り正確に）予測する ‒ 「データ→未知データ が得られる確率」を推論 汎化損失 𝐺𝑛 = − න𝑞 𝑥 log 𝑝 ∗ 𝑥|𝑋 𝑛 𝑑𝑥 を小さくする • 「構造の発見」（説明重視） ‒ データの生成源（に可能な限り近い分布）を見つける ‒ 「データセット が得られる確率」を推論 自由エネルギー 𝐹𝑛 = − log 𝑍 𝑋 𝑛 を小さくする • それぞれの定式化は上記 • これらの数式は何？？ 53

正確な予測 • 汎化損失を と定義する。 • 𝐺𝑛 = − ‫ 𝑥 𝑞 ׬‬log 𝑝∗ 𝑥|𝑋 𝑛 𝑑𝑥 = − ‫ 𝑥 𝑞 ׬‬log 𝑞 𝑥 𝑑𝑥 + 𝐾𝐿 𝑞||𝑝∗ と 変形できる • 𝑆 = − ‫ 𝑥 𝑞 ׬‬log 𝑞 𝑥 𝑑𝑥はモデルに依存しない定数 • 𝐾𝐿 𝑞||𝑝 ∗ は真の分布と予測分布のKL情報量→小さいほど予測が正確 情報理論の言葉を使うと、 真の分布と予測分布の交差（クロス）エントロピーのことです。 54

構造の発見 • データセットの真の分布：𝑄 𝑥 𝑛 = ς𝑛𝑖=1 𝑞 𝑥𝑖 ‒ 𝑥 𝑛 = 𝑥1 , … , 𝑥𝑛 ：データセット（大文字は確率変数、小文字は実現値） • データセットの確率分布を推論したい ‒ モデルから推論したこの分布が真に近いほど、 データの生成を間違いにくく説明するモデルになる • ベイズ推論では 周辺尤度𝑍 𝑋 𝑛 • やること：𝑍 𝑋 𝑛 の最大化⇔𝐹𝑛 ≔ − log 𝑍 𝑋 𝑛 の最小化 55

構造の発見 • 最大化の正当化： − න𝑄 𝑥 𝑛 log 𝑍 𝑥 𝑛 𝑑𝑥 𝑛 = − න𝑄 𝑥 𝑛 log 𝑄 𝑥 𝑛 𝑑𝑥 + 𝐾𝐿 𝑄||𝑍 ‒ 上式が成立するため、汎化損失の最小化と同様の理由で -‫ 𝑛 𝑥 𝑄 ׬‬log 𝑍 𝑥 𝑛 𝑑𝑥 𝑛 （総和損失という）の最小化が正当化される ‒ ‫𝑛𝑥 𝑄 ׬‬ 𝑑𝑥 𝑛 はデータセットの出方についての平均 → 𝐹𝑛 ≔ − log 𝑍 𝑋 𝑛 の最小化は平均的に𝐾𝐿 𝑄||𝑍 を最小化する // • 自由エネルギーを最小にするモデルはn→∞のとき真のモデルになる ‒ ただし、モデルの中に真の分布がぴったり含まれている（真のモデルと呼ぶ）必要がある ‒ たとえ真の分布がモデル候補に含まれていても、汎化損失最小化では実現できない 56

統計的推論の目的と評価量まとめ • 統計的推論の間違いは目的依存で次の量で定義する： • 目的：未知データの正確な予測 → 推論の間違い：汎化損失 • 目的：データセットの生成構造の発見・説明 → 推論の間違い：自由エネルギー • 何を「良い統計モデル」と考えるかによって i.e. 目的によって これらの評価量を使い分ける 57

FAQ • Q：この２つ以外の指標は？ • A：本来モデルの評価は目的依存で行われるべき。 最も一般的な状況では、汎化損失と自由エネルギーが重要。 • Q：真の分布qはわからないのに、指標がqに依存していない？ • A：その困難を避けるために、次ページ以降の方法を使う。 これらの方法が生まれるのに数理科学・理論研究が必要であった。 58

汎化損失と自由エネルギーの計算 • 先ほどの汎化損失は真の分布qに関する平均で定義されている ‒ 実際には真の分布qはわからないため、そのままでは計算不可能 • 真の分布は必要ないが自由エネルギーは周辺尤度という一般に は積分計算が困難な量が含まれている ‒ 自由エネルギーを計算するうえでMCMCサンプリングの平均は有効ではない 汎化損失と自由エネルギー • 𝐺𝑛 = − ‫ 𝑥 𝑞 ׬‬log 𝑝 ∗ 𝑥|𝑋 𝑛 𝑑𝑥：真の分布q(x)に依存する。 • 𝐹𝑛 = − log ‫𝑤𝑑 𝑤 𝜑 𝑤| 𝑛 𝑋 𝑃 ׬‬：積分計算が困難。 59

経験損失 • 手元のデータ𝑋 𝑛 = 𝑋1 , … , 𝑋𝑛 の損失は次式で定義される経験損失 である： • 𝑛 → ∞のとき𝑇𝑛 → 𝐺𝑛 • 経験損失で評価していい？ 60

経験損失と過学習 • 結論：経験損失で評価すると過学習を起こすのでダメ 61

経験損失と過学習 • モデルがデータを過剰に学習してしまい、未知データに対する 性能が悪化することを、over fitting（過学習または過適合）という ベイズ予測分布の過学習は 分散が大きくなる形で現れます。 62

経験損失と過学習 • 結論：経験損失で評価すると過学習を起こすのでダメ • 理由１：経験損失は手元にある「たまたま得られたデータ」についての 損失 ‒ 「見かけ上の推論の間違い」でしかない、とも言える ‒ 偶然得られたデータにぴったり当てはまっても未知データについては…… • 理由２：モデルが複雑なほど同じnのときの経験損失が小さくなる ‒ もっとも複雑なモデルが選ばれることになる ‒ 直線より放物線，放物線より三次函数…… ‒ データ点をすべて通るn次函数が選ばれることも 63

経験損失と汎化損失 • 経験損失と汎化損失は異なり，経験損失に注目すると過学習する • 汎化損失はどのような振る舞いをするのだろうか？ • 経験損失と汎化損失はどの程度異なるのだろうか？ 次の数理科学の問題に帰着例（漸近論）： １． 𝑛 → ∞のとき確率変数𝐺𝑛 の漸近挙動を求めよ． ２．確率変数𝐺𝑛 と𝑇𝑛 の差𝑉𝑛 ′もまた確率変数である． 𝑉𝑛 ′の漸近挙動を求めよ． 主要な参考文献：[Watanabe, 2009] 64

ベイズ統計理論 • 比較的明らかな結果： 𝔼 𝐺𝑛 = 𝔼 𝐹𝑛+1 − 𝐹𝑛 • 前述のとおり自由エネルギー𝐹𝑛 もまた重要な確率変数であった • まず，一般の自由エネルギーの漸近挙動が解明された 𝐹𝑛 = 𝑛𝑆𝑛 + 𝜆 log 𝑛 − 𝑚 − 1 log log 𝑛 + 𝑜𝑝 log log 𝑛 ‒ Snは経験エントロピー • 第1の主要項の係数𝜆や第2項の𝑚とは何か？ 65

ベイズ統計理論 • 一般の自由エネルギーの漸近挙動 𝐹𝑛 = 𝑛𝑆𝑛 + 𝜆 log 𝑛 − 𝑚 − 1 log log 𝑛 + 𝑜𝑝 log log 𝑛 • 係数𝜆や𝑚とは何か？ • 事後分布が正規分布で近似可能なとき（正則モデル）： ‒ この結果は古くから知られていた→ベイズ情報量規準BIC[Schwarz, 1978] • 近似不可能な場合を含む一般の場合はλは実対数閾値、mは多重度である。 66

ベイズ統計理論 • 自由エネルギーに登場したものと同じλを用いて、汎化損失 と経験損失の漸近挙動を記すことができる： 汎化損失・経験損失においても、 • 一般の汎化損失と経験損失の漸近挙動 （実現可能な） 正則モデルの場合は 𝑑 𝜆=𝜈= 2 が成立する。 こちらも古くから知られていた →赤池情報量規準AIC[Akaike, 1977] • ここで新たに登場した𝜈は特異ゆらぎと呼ばれ、対数モデルの事後分散の総和 を用いて2𝜈 = lim 𝔼 𝑉𝑛 と書ける。𝑉𝑛 は汎函数分散と呼ばれる。 𝑛→∞ 67

学習理論のゼータ函数 一般に、qからpへのKL情報量すなわち平均誤差函数をK、事前分布をφとする とき、学習理論のゼータ函数が次で定義される： ζ(z)=∫K(w)zφ(w)dw. • • ζ(z)は複素数平面全体に有理型函数として一意に解析接続され、その極 はすべて負の有理数となることが証明できる[Atiyah, 1970]。 最大極の絶対値を実対数閾値と，その極の位数を多重度と呼ぶ。 平均誤差関数は𝑲 𝒘 = ‫𝐠𝐨𝐥 )𝒙(𝒒 ׬‬ 𝒒𝒙 𝒑 𝒙𝒘 𝒅𝒙 具体的なモデルについてFn及びGnの挙動を知りたければ、 上記ゼータ函数の最大極を調べればよい。 →たくさんの研究(Aoyagi, Drton, Rusakov, Yamazaki,…) 𝟏/ 𝒛 + 𝝀 𝐗 𝐗 𝐗𝐗 𝑚 𝐗 𝐎 ℂ 公式10本ノック？watanabe-www.math.dis.titech.ac.jp/users/swatanab/joho-gakushu5.html 68

• http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/joho-gakushu5.html

ベイズ統計理論の応用 • 上記のように、汎化損失や自由エネルギーの漸近挙動が得られた • モデルの評価にはどう使えばいいか？ ‒ 学習理論のゼータ函数にも真の分布q(x)は登場してしまう……。 ‒ 真の分布の仮定ごとにゼータ函数を解析？ 厳しい……。 ‒ 漸近論は役に立たない？ • 漸近挙動を用いて、新しい情報量規準が作られた。 69

情報量規準 • データから簡便に計算できる量を使って汎化損失や自由エネルギーを計算 できる量はしばしば情報量規準と呼ばれる。 1 𝑑 1 𝑛 𝑛 𝑛 ‒ 赤池情報量規準:AIC = − log 𝑃 𝑋 𝑛 |ෝ 𝑤𝑀𝐿𝐸 + , 𝔼 AIC = 𝔼 𝐺𝑛 + 𝑜 𝑑 ‒ ベイズ情報量規準:BIC = − log 𝑃 𝑋 𝑛 |ෝ 𝑤𝑀𝐿𝐸 + log 𝑛 , BIC = 𝐹𝑛 + 𝑂𝑝 1 2 • AICやBICは正則モデルのみで妥当性を持つ。一般の場合は？ →前述の漸近理論から新たな情報量規準が作られた。 ‒ ベイズ推論の漸近理論を用いているため、ベイズ推論において利用可能 70

広く使える情報量規準 ベイズ統計理論の結果から、新たな情報量規準が作られた。 • 広く使える情報量規準WAIC[Watanabe, 2010]: • 経験損失𝑇𝑛 と汎函数分散𝑉𝑛 は共にデータから計算できる量である。 71

広く使える情報量規準 前述の漸近挙動より、 𝜆 𝜆−2𝜈 𝑛 𝑛 𝔼 𝐺𝑛 − 𝔼 𝑇𝑛 = 𝑆 + − 𝑆 − ‒ 実は𝑜 1 𝑛 は𝑜 1 𝑛2 +𝑜 1 𝑛 = 2𝜈 𝑛 +𝑜 1 𝑛 = 𝔼 𝑉𝑛 𝑛 +𝑜 1 𝑛 にまで改善可能であり、 WAICは1個抜き交差検証（LOOCV）損失と漸近等価である。 72 .

WAICの計算 WAICは次のように計算することができる。 1. 事後分布𝜓 𝑤 𝑋 𝑛 からサンプリングする。 log 𝑝 𝑋𝑖 𝑤𝑘 を (i ,k) 成分とする テーブルがあれば計算可能！ 𝑤𝑘 ∼ 𝜓 𝑤 𝑋 𝑛 , 𝑘 = 1, … , 𝐾 2. 経験損失𝑇𝑛 と汎函数分散𝑉𝑛 を計算する。 1 1 n 𝐾 𝑇𝑛 ← − σ𝑛𝑖=1 log σ𝐾 𝑘=1 𝑝 𝑋𝑖 𝑤𝑘 𝑉𝑛 ← 3. 1 σ𝑛𝑖=1 σ𝐾 𝐾 𝑘=1 log 𝑝 𝑋𝑖 𝑤𝑘 , 2 − 1 σ𝐾 𝐾 𝑘=1 log 𝑝 𝑋𝑖 𝑤𝑘 2 ２．で計算した量からWAICが得られる。 WAIC ← 𝑇𝑛 + 𝑉𝑛 𝑛 73

広く使えるベイズ情報量規準 ベイズ統計理論の結果から、新たな情報量規準が作られた。 • 広く使えるベイズ情報量規準WBIC[Watanabe, 2013]: 𝛽 ‒ 𝔼𝑤 ⋅ は「逆温度βの事後分布𝜓𝛽 」による平均: • 対数尤度log 𝑃 𝑋 𝑛|𝑤 も𝜓𝛽 からのサンプリングもデータから計算できる。 74

広く使えるベイズ情報量規準 前述の漸近挙動より、𝐹𝑛 = 𝑛𝑆𝑛 + 𝜆 log 𝑛 − 𝑚 − 1 log log 𝑛 + 𝑂𝑝 1 . 𝛽 一方、β0=1/log(n)のとき−𝔼𝑤0 log 𝑃 𝑋 𝑛|𝑤 = 𝑛𝑆𝑛 + 𝜆 log 𝑛 + 𝑂𝑝 𝜆 log 𝑛 であることを証明できる。 ‒ 自由エネルギーは右図の青曲線の下面積であり、 自由エネルギー モデルに依存せず 1/log(n) のときの対数尤度の 負の事後平均（WBIC）に漸近的に等しい。 WBIC ‒ 漸近挙動の主要項まで一致し、データから計算しやすい 量としてWBICが与えられる。 β=0 β=1 β=1/log(n) 75 β

WBICの計算 WBICは次のように計算することができる。 1. 逆温度β=1/log(n)の事後分布𝜓𝛽 𝑤 𝑋 𝑛 からサンプリングする。 𝛽 𝑤𝑘 ∼ 𝜓𝛽 𝑤 𝑋 𝑛 , 𝑘 = 1, … , 𝐾 ‒ メトロポリスヘイスティング法なら𝐻 𝑤 = − log 𝑃 𝑋 𝑛 𝑤 − log 𝜑 𝑤 を 𝐻𝛽 𝑤 = −𝛽 log 𝑃 𝑋 𝑛 𝑤 − log 𝜑 𝑤 に置き換えればよい。 2. 分布𝜓𝛽 について対数尤度を平均し、符号反転させることでWBICが得られる。 WBIC = 1 − σ𝐾 log 𝑃 𝐾 𝑘=1 𝑋 𝑛 𝛽 |𝑤𝑘 1 𝛽 𝑛 σ = − σ𝐾 𝑘=1 𝑖=1 log 𝑝 𝑋𝑖 𝑤𝑘 𝐾 𝛽 log 𝑝 𝑋𝑖 𝑤𝑘 を (i ,k) 成分とする テーブルがあれば計算可能！ 76

情報量規準の使い方 統計モデリングで直面する状況は例えば以下の状況である。 • 検討するモデルと事前分布の候補が複数ある。その中から最も妥当 なものを選択したい。 1. 未知データに対する予測を最も間違わずに行えるモデルはどのモデルか？ 2. データを発生している分布がモデルに含まれるとしたらどのモデルか？ 状況１では汎化損失を、状況２では自由エネルギーを最小化したい。 しかし真の分布は不明であり、自由エネルギーを直接計算することは 難しい。そこで、状況１ではWAICを、状況２ではWBICを計算する。 77

情報量規準の使い方 情報量規準を用いたモデル選択は以下の流れとなる。 1. 検討対象の複数のモデルそれぞれについてベイズ推論を行う。 ‒ 事後分布・逆温度βの事後分布・予測分布のうち必要なものを計算する。 2. それぞれの推論結果に対してWAICやWBICを計算する。 3. 目的に応じて、WAICやWBICを最小化するモデルを選択する。 ‒ 未知データへの予測⇔汎化損失最小化”⇔”WAIC最小化 ‒ データ発生源の発見⇔自由エネルギー最小化”⇔”WBIC最小化 WAICは平均値と分散が汎化損失のそれと 漸近等価ですが逆相関です。 WBICは主要項のみが一致する等価性で、 数値的な分散が小さくありません。 “⇔”はこの意味での近似的な⇔です。 78

まとめ 本チュートリアルではベイズ統計学について以下を扱った。 • ベイズ推論の基礎概念と推論フレームワーク • 代表的な推論アルゴリズム（ﾒﾄﾛﾎﾟﾘｽﾍｲｽﾃｨﾝｸﾞ、ｷﾞﾌﾞｽｻﾝﾌﾟﾘﾝｸﾞ） • ベイズ推論の汎化・経験損失や自由エネルギーの一般理論と、 それを用いたモデル選択指標である情報量規準WAIC・WBIC 79

参考文献 [Akaike, 1977]: “Likelihood and Bayes Procedure” [Atiyah, 1970]: “Resolution of singularities and division of distributions” [Schwarz, 1978]: “Estimating the dimension of a model” [Watanabe, 2009]: “Algebraic Geometry and Statistical Learning Theory” [Watanabe, 2010]: “Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory” [Watanabe, 2013]: “A widely applicable Bayesian information criterion” 80