カイヤン雑記帳

カイヤンがやったことを書いておいたり、ぼやきたいことを書き込んだりする場所

ベイズ統計学の逆温度パラメータ

おはようございますまたはこんにちはまたはこんばんは,カイヤンです.

就活エントリーが有名な方に捕捉されたりものすごくバズったりと過分な評価をいただいており,とても驚いているこの頃です.

私や私のドクトリンが評価されているのではなく,記事から就活最前線の様子を推測しやすいからこれだけ評価をいただいているに過ぎないと肝に銘じて思い上がらぬよう過ごしたいです *1

さて,今回はPeingに届いた質問に関して短く書いてみようと思います.

ベイズ推定における逆温度パラメータって具体的に何を制御しているものですか?(最尤法とのつながりなどでよく出てきま | Peing -質問箱-

という質問をいただきました.

質問箱では簡単に答えていますが,もう少し語らせてください.

筆者質問箱の回答

「対数尤度+正則化項」という損失を考えると結局正則化の効果度(逆温度が大きいほど弱い正則化,無限大への極限で最尤法)と捉えることになりますが,そもそもベイズ推定だと損失最小化推定量が目的でないので「事後分布の形状」を司るということですかね

これではハイコンテクスト過ぎるので,もう少しかみ砕いていきます.

ベイズ推測の逆温度

そもそも質問にあったベイズ推測の逆温度ですが,統計学機械学習を学んでいても聞きなれない言葉かもしれませんので,まずその定義から説明します.

ベイズ推測を行うにあたって,確率モデル p(x|w)と事前分布 \varphi(w)を用意します.ここで xはデータの空間の元で wはパラメータの空間のそれです. このとき n個のデータ  X^n = (X_1,\ldots,X_n)が得られているときの事後分布は

 \psi(w|X^n):=\frac{\prod_{i=1}^n p(X_i|w) \varphi(w)}{Z_n(X^n)}

で定義されるのでした.ここで,

 P(X^n|w) = \prod_{i=1}^n p(X_i|w)

を尤度といい,

 Z_n(X^n) = \int P(X^n|w)\varphi(w) dw

を周辺尤度と言います.この事後分布を用いて予測分布

 p^*(x) = \int p(x|w)\psi(w|X^n) dw

が定義されますが,この「予測分布がデータを発生している分布であろう」と推測することがベイズ推測です.

これに対して,逆温度 \betaの事後分布は次の式で定義されます.

 \psi(w|X^n,\beta) := \frac{\prod_{i=1}^n p(X_i|w)^{\beta} \varphi(w)}{Z_n(X^n,\beta)}

ここでも Z_nは事後分布の積分が1になるようにするための正規化定数です.

逆温度,というと物理的な意味を考えたくなるかもしれません.実際その通りでこれは完全に統計力学からの輸入です. 私は統計力学については全くの素人なのですが,「ハミルトニアン H(w)が確率過程であるすなわちランダムハミルトニアンを考えているときのボルツマン分布

 p(w) = \frac{1}{Z(\beta)}\exp(-\beta H(w))

が逆温度 \betaの事後分布であり,分配函数 Z(\beta)について Z(1)が周辺尤度である」という対応があるようです. H(w)(の実現値)はエネルギーなので逆温度(ボルツマン定数と温度の積の逆数→エネルギーの逆数と同じ次元量)をかけて無次元化して指数部分となっています(確率は無次元量です).

統計学の範囲では,「上記の逆温度 \betaの事後分布に出てくる制御変数 \betaを逆温度と呼ぶ」,ということで問題ないはずです. 標語的には,「統計学では逆温度がぴったり1であるときが常温」と言うこともできるでしょう*2

制御変数「逆温度」を変えると何が起こるか

前節でベイズ推測の逆温度は事後分布の制御変数であると言いました.「事後分布を制御している」というのではあまりにもざっくりしすぎているため*3,掘り下げていきましょう.

最尤推定・事後確率最大化推定と逆温度

まず,質問箱にもあったように最尤推定と絡めて見てみます.

最尤推定では,「パラメータはある定数であり,尤度を最大化するパラメータがその求めるべき”ある定数”であろう」という考え方で行われるものです. 数式の上では上述の Pを最大化しますが確率のサンプルサイズ個の積という小さくなりがちな値を扱う関係で負の対数を取った

 L_n(w) := -\log P(X^n|w) = - \sum_{i=1}^n \log p(X_i|w)

の最小化を行うことが多いです.

 L_n(w)の最小化とはすなわち実函数の最小化ですから,高校数学で習った微分や大学1年生で扱った多変数関数の微分が役に立ちます. そのときの試験問題のように解析的に最小化元を求められないことは実際のモデルではしばしば起こるので,微分を計算して勾配降下法などにより数値的に求めるということもよく行われるものです. 負の対数尤度の最小化元 \hat{w}_{MLE}のことを最尤推定量と言います.

一方で,事後確率最大化推定は「パラメータは確率変数であり,事前分布をデータに基づいて更新した事後分布がその分布であろう」という点でベイズ推測と同じです. しかし,ベイズ推測では予測分布を以てデータを発生している真の分布を推測するのに対して,事後確率最大化推定では「得られた事後分布を最大化するパラメータを点推定した結果を確率モデルに戻したものが真の分布であろう」と考えます.分布を推測するベイズ推測とは異なり,推定量を求めているという点で最尤推定に近いです.

日本語で色々と書いていますが,結局のところ事後確率最大化推定は負の対数事後分布

 L_n(w) := -\log P(X^n|w) \varphi(w) = - \sum_{i=1}^n \log p(X_i|w) - \log \varphi(w)

の最小化元を求めるということです*4. これも最適化問題として最尤推定と同じように解かれます.得られる推定量は,事後確率最大化推定量と呼ばれます.

このように,推定量の計算にあたって最小化される函数はしばしば損失函数と呼ばれます. 損失函数を最小化する推定量を求める(点推定する,とも言います)という意味で,最尤も事後確率最大化も等価です. 実際,上の2つの損失函数を見てみると,事前分布に関する項があるかないかという違いしかありません. ここで例えば事前分布として平均0の正規分布を用いると事前分布の項はリッジ正則化項に,平均0のラプラス分布を用いるとLASSO正則化項にそれぞれなります. この場合では事前分布の分散の逆数が正則化パラメータ lに対応します:

 L_n(w) =- \sum_{i=1}^n \log p(X_i|w) - l \|w\|^2.

さて,事後確率最大化推定を逆温度が1とは限らない事後分布について行ってみましょう.最小化すべき損失函数

 L_n(w) := - \sum_{i=1}^n \beta\log p(X_i|w) - \log \varphi(w)

となりますが, wに依存しない定数倍を施しても最小化元は変わらないため 1/\beta倍した

 L_n(w) := - \sum_{i=1}^n \log p(X_i|w) - \frac{1}{\beta}\log \varphi(w)

の最小化と等価です.

事前分布の項に逆温度の逆数がかかっています.ここで \beta=1とすると元の事後確率最大化推定になります. 一方で, \betaを大きくしていくと正則化という観点では弱い正則化になります.逆に小さな逆温度では正則化は強くなります. つまり,点推定を考えているときは,逆温度は正則化の強弱を制御していることがわかります.

注意

よく上の事後確率最大化の損失函数を見て,「事前分布が一様分布であるとき,事後確率最大化推定と最尤推定が一致する」と言われます. これは最小化する問題としては等価ではありますが,そもそも最尤推定では「パラメータはある一点であって揺らいでいない」ということになるため厳密には 「事前分布が最尤推定量に質量をもつデルタ分布 \delta(w-\hat{w}_{MLE})であるとき,事後確率最大化推定と最尤推定が一致する」ということになります.

(2019/5/8追記)$\beta \rightarrow \infty$としたとき事後分布がどうなるかを調べている方がいました:

t.co

最尤推定は主観依存の事前分布を用いないのでベイズ推測や事後確率最大化推定よりも正しい!」という意見が未だに出てきますが, 実際は最尤推定では「最尤推定量に質量をもつデルタ分布」という極めて極端な事前分布を使っており,ベイズと同程度に客観的ではありません. なお,統計学は不良設定問題*5を扱う学問であるため,そもそも「正しい統計的推測」というものはこの世に存在しないことは 20世紀後半からは統計学を学んだ人にとってはよく知られていることと思います.

ベイズと最尤のどちらが”正しい”のか」「事前分布は何が正しいのか」「統計学は不良設定問題を扱う学問である」「正しい統計的推測は存在しない」といった事柄に関心を持たれた読者の皆様は以下のホームページ・講演資料を参考にされると幸せになれると思います.

渡辺澄夫,ベイズ推論:いつも何度でも尋ねられること

渡辺澄夫,電子情報通信学会ソサイエティ大会, AI-2 データ科学とコンピュータ科学の基礎理論と展開, 2016年9月20日北海道大学

ベイズ推測

前小節では,パラメータを点推定する場合で逆温度が何を制御するかを見ていました.それでは予測分布を推測するベイズ推測ではどうでしょうか.

結論を先に書いてしまうと,事後分布の形状を制御しています.ここで形状というのは,分布の峰の鋭さや漸近的な性質などです.

一般に事後分布は正規分布などの簡単な分布で近似することはできませんので,峰といっても複数あることは少なくありません. そのような複雑な確率分布を解析的に計算することは困難です.そこで,実際にベイズ推測を行うときは本記事冒頭で述べた事後分布は近似的に求めることが多いです. いくつか近似推測がありますが,その中でもマルコフ連鎖モンテカルロ法(MCMC)と呼ばれる方法がメジャーでしょう.この方法を使うと理論上は無限にサンプリングをすることで真の事後分布を実現できるためです. 実際には無限にサンプリングすることはできないため近似となります.

さて,サンプリングと言いますが,事後分布からまんべんなくサンプリングできないと意味がありません. しかしサンプリング手法によってはそれは意外と難しいことです.パラメータをいくつかに分けてそれぞれの分布を簡単に表すことができると, ギブスサンプリングを使えますが,これは使えるモデルや事前分布が限られます.一般には,事後分布をエネルギー函数――上で述べたランダムハミルトニアンです――に沿ってサンプリング箇所が動き回るようにします. 例えばメトロポリスヘイスティング法はそのためのプリミティブかつ有効な手法です.しかし,事後分布そのものはすぐ上で言ったように多峰性を持っていることが多く, これがとても深い場合は谷底にはまってしまうという問題が起きます.

(注)MCMCは確かに多変数関数のある種の勾配に沿って移動しますが,何かを探しているわけではなくそこらじゅうを適度に移動させることが目的です.そのため,MCMCが収束するというのは一つの確率分布に従うサンプル列に収束するということです.この確率分布を定常分布といい,定常分布が実現したい事後分布になるようにMCMCを設計します.

サンプリング箇所が局在してしまうと,MCMCによるサンプルは偏ってしまい事後分布からのサンプリングとみなせなくなります. このことは深い谷底にはまってしまうことが所以です.メトロポリスヘイスティング法では, t番目のサンプリング箇所 w_tから次の箇所 w'へ移動するかどうかを確率により決めます. その確率は

 a(t) = \min \{ 1, \exp(-(H(w_t) - H(w')))

で定義され,遷移確率と言います.確率 a(t)で次のサンプリング箇所 w_{t+1} w'に,確率 1-a(t)でその場に留まって w_tとなります. 候補と現在位置のエネルギー差が大きい場合は遷移確率は0に近くなってしまうため,まんべんなくサンプリングすることができません.

この問題を回避するために,逆温度が使われます:逆温度を変えることで事後分布の形状を変化させ,谷底を浅くするのです. それは,レプリカ交換モンテカルロ法(レプリカ法)と呼ばれる手法です.一般に,逆温度 \betaの事後分布を考え,様々な逆温度を利用して逆温度1の事後分布からのサンプリングを効率的に行います.

 \psi(w|X^n,\beta)=\frac{1}{Z(\beta)}\exp(-\beta H(w))

という事後分布で,逆温度は0以上1以下のいくつかの値とします.このいくつかの逆温度それぞれで事後分布を考え, 普通にメトロポリスヘイスティング法などでサンプリングを行います.そしてサンプリングの最中にある確率で逆温度を交換します. 例えば,逆温度 \beta_1でサンプリングして

w^1_1,\ldots,w^1_L

というサンプルを,逆温度 \beta_2でサンプリングして

w^2_1,\ldots,w^2_L

というサンプルをそれぞれ得ていた時に ある確率で次の L個を合わせたサンプルが

 w^1_1,\ldots,w^1_L,w^2_{L+1},\ldots,w^2_{2L}

及び

 w^2_1,\ldots,w^2_L,w^1_{L+1},\ldots,w^1_{2L}

になるということです.各 L個のサンプルがそれぞれレプリカと呼ばれ,それが確率的に交換されながらサンプリングが進んでいきます. まんべんなくサンプリングするためには,各温度で交換確率が一定であることが重要ですが, そのためには逆温度の刻みを等比数列になるように設計すればよいことが知られています. 詳細は

  • 永田賢二,"Theory and Experiments of Exchange Ratio for Exchange Monte Carlo Method", 2008

をご覧ください. このレプリカの交換確率についても面白い数理が隠れている*6のですが,それも上記文献にあります.結果だけであれば,

  • 渡辺澄夫,『ベイズ統計の理論と方法』,2012

にもあります.

むすび

以上のように,逆温度は事後分布の形状を制御しており,事後分布の近似推論に活用されています. 事後分布の形状という点から再び事後確率最大化推定を考えると,逆温度が無限大の極限では事後分布は最尤推定量に質量をもつデルタ分布になっていると考えられます.

短めといいつつ思ったより長くなってしまいました(レプリカ交換モンテカルロは難しいのう……). ここまで読んでいただき,ありがとうございました! またどこかの記事でお会いできれば光栄です.

*1:このブログの閲覧数なんて1日50も行けば多い方なのに落ち着き気味の週明け29日ですら300超でピークは2000弱とか興奮せざるを得なかった.

*2:分野により「常温」が異なるというネタの統計学版です(笑).ボルツマン定数のことを考えると,これはこのページをご覧になっている皆様がおいでのお部屋などと比べとんでもない高温です!

*3:統計モデルによっては事前分布のパラメータ=ハイパーパラメータによって事後分布が大きく異なることがあります.

*4:周辺尤度はパラメータに依存しないため省略できます.むしろ周辺尤度は解析的に計算できないことが多いため出てこられると困ります.

*5:基礎となる確率すなわちデータを発生している分布が未知

*6:ここにも代数幾何学との交わりが出てきます:実対数閾値が厳密な交換確率を与えます!