自由エネルギーにかんするノート

《自由エネルギー》は統計モデルを使ったデータで，モデル評価の指標として使われる推定量の一つです．実用上はStanなどを使って，MCMCサンプルからこれらの指標を計算することが多いと思います．

その際，解析的に計算した厳密値とPCで計算した値が一致するかどうか気になったことはないでしょうか？

そこで本ノートでは，《解析的に計算した（手計算で求めた）自由エネルギー》と《PCを使って計算した近似値》が一致するのかどうかを単純な例を使って調べます ^1 ．以下，数式の例は渡辺（2012）を，計算コードの例は浜田・石田・清水（2019）を参照しました．

自由エネルギー

まずは定義を確認します．

確率モデルを$p(x|w)$，事前分布を$\phi (w)$とおき，事後分布を
$$p(w|x^n)=\frac{1}{Z_n}\prod _{i=1}^{n}p(x_i|w)\phi (w)$$と定義します．ここで関数$Z_n$は次の積分
$$Z_n=\int \prod _{i=1}^{n}p(x_i|w)\phi (w)dw$$として定義され，分配関数や周辺尤度(marginal likelihood)と呼ばれます．

ここで，周辺尤度の対数の符号を逆転させたものを，自由エネルギー(free energy)と定義し、記号$F_n$で表します.
$$\begin{array}{rcl}& F_n& =-\log Z_n\\ & & =-\log \int \prod _{i=1}^{n}p(x_i|w)\phi (w)dw\end{array}$$
周辺尤度$Z_n$は，サンプル$X^n$についての分布となっており（$x^n$で積分すると1になる），私たちがモデルとして設定した確率分布とみなすことができます．

指数分布族の分配関数・自由エネルギー・事後分布

一般的に確率モデルは
$$p(x|w)=v(x)\exp (f(w)\cdot g(x)),f(w),g(x)\in {\mathbb{R}}^k$$という形で書けるとき指数分布族と言います（$f(w)\cdot g(x)$はベクトル値関数の内積）．

ここで$w$の事前分布として
$$\phi (w|\varphi )=\frac{1}{z(\varphi )}\exp (\varphi \cdot f(w))$$という形の事前分布を考えます（$\varphi \in {\mathbb{R}}^k$をハイパーパラメータなどと呼びます）．この事前分布を《共役な事前分布》と言います．

事前分布を$w$で積分すると1になるので
$$z(\varphi )=\int \exp (\varphi \cdot f(w))dw$$ です．

共役事前分布を使って，確率モデルが指数分布族である場合の分配関数を一般的に計算します．
$$\begin{array}{rcl}& Z_n& =\int \prod _{i=1}^{n}p(x_i|w)\phi (w)dw\\ & & =\int \prod _{i=1}^{n}p(x_i|w)\frac{1}{z(\varphi )}\exp (\varphi \cdot f(w))dw\\ & & =\frac{1}{z(\varphi )}\int \prod _{i=1}^n\{v(x_i)\exp (f(w)\cdot g(x_i))\}\exp (\varphi \cdot f(w))dw\\ & & =\frac{1}{z(\varphi )}\int (\prod _{i=1}^{n}v(x_i))\prod _{i=1}^n\exp (f(w)\cdot g(x_i))\exp (\varphi \cdot f(w))dw\\ & & =\frac{1}{z(\varphi )}\int (\prod _{i=1}^{n}v(x_i))\prod _{i=1}^n\exp (g(x_i)\cdot f(w)+\varphi \cdot f(w))dw\\ & & =\frac{1}{z(\varphi )}\int (\prod _{i=1}^{n}v(x_i))\prod _{i=1}^n\exp ((g(x_i)+\varphi )\cdot f(w))dw内積の性質より\\ & & =\frac{1}{z(\varphi )}\int (\prod _{i=1}^{n}v(x_i))\exp \{(\varphi +\sum _{i=1}^{n}g(x_i))\cdot f(w)\}dw\end{array}$$

ここで共役事前分布の仮定より $$\hat{\varphi }=\varphi +\sum _{i=1}^{n}g(x_i)$$とおけば

$$\begin{array}{rclr}& z(\hat{\varphi })& =\int \exp (\hat{\varphi }\cdot f(w))dw& \\ & z(\varphi +\sum _{i=1}^{n}g(x_i))& =\int \exp ((\varphi +\sum _{i=1}^{n}g(x_i))\cdot f(w))dw\end{array}$$です．

これを使えば分配関数は
$$Z_n=(\prod _{i=1}^{n}v(x_i))\frac{z(\hat{\varphi })}{z(\varphi )}$$
です．ただし$$\hat{\varphi }=\varphi +\sum _{i=1}^{n}g(x_i)$$です．

自由エネルギー$F_n$は分配関数$Z_n$の対数の符号を反転した量だから
$$\begin{array}{rc}{F}_n& =-\log \left((\prod _{i=1}^{n}v(x_i))\frac{z(\hat{\varphi })}{z(\varphi )}\right)\\ & =-\sum _{i=1}^n\log v(x_i)-\log \frac{z(\hat{\varphi })}{z(\varphi )}\end{array}$$です． また事後分布は
$$\begin{array}{rcl}& p(w|x^n)& =\frac{1}{Z_n}\prod _{i=1}^{n}p(x_i|w)\phi (w)\\ & & =\frac{1}{Z_n}\frac{1}{z(\varphi )}(\prod _{i=1}^{n}v(x_i))\exp (\hat{\varphi }\cdot f(w))\\ & & ={\left((\prod _{i=1}^{n}v(x_i))\frac{z(\hat{\varphi })}{z(\varphi )}\right)}^{-1}\frac{1}{z(\varphi )}(\prod _{i=1}^{n}v(x_i))\exp (\hat{\varphi }\cdot f(w))\\ & & =\frac{1}{z(\hat{\varphi })}\exp (\hat{\varphi }\cdot f(w))\\ & & =\phi (w|\hat{\varphi })=\phi (w|\varphi +\sum _{i=1}^{n}g(x_i))\end{array}$$です．

この計算から，《指数分布族に対しては共役事前分布を定義でき，事後分布は共役事前分布と同じ関数となること》が分かります．ただし事前分布のハイパーパラメータは$\varphi$から$\hat{\varphi }=\varphi +\sum _{i=1}^{n}g(x_i)$に変化します．$g(x_i)$は確率モデルのカーネルですから，事後分布（のハイパーパラメータ）は，データを確率モデルというフィルターを介して学習した結果と見なせるでしょう．

正規分布の共役事前分布

正規分布のパラメータ$\mu$の共役事前分布を作ってみましょう．正規分布の確率モデル$p(x|w)$は
$$\begin{array}{rcl}& p(x|w)& =\mathrm{Normal}(x|\mu ,{\sigma }^2)\\ & & =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \{-\frac{(x-\mu {)}^2}{2{\sigma }^2}\}\\ & & =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \{-\frac{x^2}{2{\sigma }^2}\}\exp \{\frac{\mu x}{{\sigma }^2}-\frac{{\mu }^2}{2{\sigma }^2}\}\\ & & =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \{-\frac{x^2}{2{\sigma }^2}\}\exp \{(\mu /{\sigma }^2,-{\mu }^2/2{\sigma }^2)\cdot (x,1)\}\end{array}$$と書けるので，$$v(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \{-\frac{x^2}{2{\sigma }^2}\}，f(w)=(\mu /{\sigma }^2,-{\mu }^2/2{\sigma }^2),g(x)=(x,1)$$ とみなすと
$$p(x|w)=v(x)\exp (f(w)\cdot g(x))$$ なので，確かに指数分布族です．

指数分布族の共役事前分布は
$$\phi (w|\varphi )=\frac{1}{z(\varphi )}\exp (\varphi \cdot f(w))$$という形で定義可能でした．
$$\varphi =({\varphi }_1,{\varphi }_2),f(w)=(\mu /{\sigma }^2,-{\mu }^2/2{\sigma }^2)$$を代入すると $$\begin{array}{rc}\phi (w|\varphi )& =\frac{1}{z(\varphi )}\exp (({\varphi }_1,{\varphi }_2)\cdot (\mu /{\sigma }^2,-{\mu }^2/2{\sigma }^2))\\ & =\frac{1}{z(\varphi )}\exp (\frac{\mu {\varphi }_1}{{\sigma }^2}-\frac{{\mu }^2{\varphi }_2}{2{\sigma }^2})\end{array}$$です．ここで$w=(\mu ,{\sigma }^2)$の${\sigma }^2$は定数であると見なし，$\phi (w|\varphi )$を$\mu$だけの分布と考えれば，$$\phi (\mu |\varphi )=\frac{1}{z(\varphi )}\exp (\frac{\mu {\varphi }_1}{{\sigma }^2}-\frac{{\mu }^2{\varphi }_2}{2{\sigma }^2})$$です．$\phi (\mu |\varphi )$が$\mu$の分布であるためには積分して1になる必要があるので
$$\int \phi (\mu |\varphi )d\mu =1$$ より
$$z(\varphi )=\int \exp (\frac{\mu {\varphi }_1}{{\sigma }^2}-\frac{{\mu }^2{\varphi }_2}{2{\sigma }^2})d\mu$$です．
ガウス積分を使って計算すると $$\begin{array}{rcl}& z(\varphi )& =\int \exp (\frac{\mu {\varphi }_1}{{\sigma }^2}-\frac{{\mu }^2{\varphi }_2}{2{\sigma }^2})d\mu \\ & & =\int \exp (-\frac{{\varphi }_2}{2{\sigma }^2}{(\mu -\frac{{\varphi }_1}{{\varphi }_2})}^2+\frac{{\varphi }_1^2}{2{\varphi }_2{\sigma }^2})d\mu \\ & & =\exp \left(\frac{{\varphi }_1^2}{2{\varphi }_2{\sigma }^2}\right)\int \exp (-\frac{{\varphi }_2}{2{\sigma }^2}{(\mu -\frac{{\varphi }_1}{{\varphi }_2})}^2)d\mu \\ & & =\exp \left(\frac{{\varphi }_1^2}{2{\varphi }_2{\sigma }^2}\right)\sqrt{2\pi {\sigma }^2/{\varphi }_2}\end{array}$$したがって$\mu$の共役な事前分布は
$$\phi (\mu |{\varphi }_1,{\varphi }_2)={(\exp \left(\frac{{\varphi }_1^2}{2{\varphi }_2{\sigma }^2}\right)\sqrt{2\pi {\sigma }^2/{\varphi }_2})}^{-1}\exp (\frac{\mu {\varphi }_1}{{\sigma }^2}-\frac{{\mu }^2{\varphi }_2}{2{\sigma }^2})$$です． 一見複雑に見えますが${\varphi }_1=0,{\varphi }_2=1$を仮定すれば，
$$\begin{array}{rcl}& \phi (\mu |{\varphi }_1,{\varphi }_2)& =\phi (\mu |0,1)\\ & & ={(\exp \left(\frac{0^2}{2{\sigma }^2}\right)\sqrt{2\pi {\sigma }^2})}^{-1}\exp (\frac{0}{{\sigma }^2}-\frac{{\mu }^2}{2{\sigma }^2})\\ & & =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp (-\frac{{\mu }^2}{2{\sigma }^2})\end{array}$$
なので，$\mu$の共役な事前分布は正規分布であることが分かります．

自由エネルギーの具体例

具体的な確率モデルと事前分布を使って，自由エネルギーを計算してみましょう．以下の例は渡辺（2012）の例をさらに簡略化したものです．

計算のための具体例． $$\begin{array}{rc}確率モデル p(x|\mu ,1)& =\mathrm{Normal}(x|\mu ,1)\\ 事前分布 \phi (\mu )& =\mathrm{Normal}(\mu |0,1)\end{array}$$確率モデルが平均$\mu$，分散1の正規分布で，事前分布が$\mu$が平均0，分散1の正規分布に従うと仮定します．

さきほどの計算結果から，この事前分布は共役事前分布であることが分かります．

確率モデルが指数分布族で，事前分布が共役事前分布であるとき，自由エネルギーは
$$F_n=-\sum _{i=1}^n\log v(x_i)-\log \frac{z(\hat{\varphi })}{z(\varphi )}$$と書けるのでした． 仮定より，確率モデルは
$$\mathrm{Normal}(x|\mu ,1)=\frac{1}{\sqrt{2\pi }}\exp \{-\frac{(x-\mu {)}^2}{2}\}$$なので，$$v(x_i)=\frac{1}{\sqrt{2\pi }}\exp \{-\frac{x_i^2}{2}\}$$です．
$$z(\varphi )=\exp \left(\frac{{\varphi }_1^2}{2{\varphi }_2{\sigma }^2}\right)\sqrt{2\pi {\sigma }^2/{\varphi }_2}$$なので，${\varphi }_1=0,{\varphi }_2=1$および${\sigma }^2=1$を代入すると
$$z(\varphi )=\sqrt{2\pi {\sigma }^2}=\sqrt{2\pi }$$ です． また$\hat{\varphi }=\varphi +\sum _{i=1}^{n}g(x_i)$と定義していたのでベクトルで書けば，
$$\begin{array}{rcl}& \hat{\varphi }& =\varphi +\sum _{i=1}^{n}g(x_i)\\ & (\widehat{{\varphi }_1},\widehat{{\varphi }_2})& =({\varphi }_1,{\varphi }_2)+(\sum _{i=1}^n{x}_i,\sum _{i=1}^{n}1)\end{array}$$なので
$$\widehat{{\varphi }_1}={\varphi }_1+\sum _{i=1}^n{x}_i,\widehat{{\varphi }_2}={\varphi }_2+\sum _{i=1}^{n}1={\varphi }_2+n$$となり $$\begin{array}{rcl}& z(\hat{\varphi })& =\exp \left(\frac{({\varphi }_1+\sum _{i=1}^n{x}_i{)}^2}{2({\varphi }_2+n)}\right)\sqrt{2\pi /({\varphi }_2+n)}\\ & & =\exp \left(\frac{(\sum _{i=1}^n{x}_i{)}^2}{2(1+n)}\right)\sqrt{2\pi /(n+1)}\end{array}$$です．
ゆえに仮定のもとで，自由エネルギーは $$\begin{array}{rcl}& F_n& =-\sum _{i=1}^n\log v(x_i)-\log \frac{z(\hat{\varphi })}{z(\varphi )}\\ & & =-\sum _{i=1}^n\log (\frac{1}{\sqrt{2\pi }}\exp \{-\frac{x_i^2}{2}\})-\log \left(\frac{\exp \left(\frac{(\sum _{i=1}^n{x}_i{)}^2}{2(1+n)}\right)\sqrt{2\pi /(n+1)}}{\sqrt{2\pi }}\right)\\ & & =-\sum _{i=1}^n(-\frac{x_i^2}{2}-\log \sqrt{2\pi })-\exp \left(\frac{(\sum _{i=1}^n{x}_i{)}^2}{2(1+n)}\right)+\log \sqrt{n+1}\end{array}$$です．

かなりシンプルな仮定から出発したおかげで，自由エネルギーを陽表的な関数として示すことができました．ここから，モデルが複雑になると，解析的な計算が困難であることは容易に想像できるでしょう．そのため実用上はMCMCサンプルを用いた近似計算がよく使われます．

パラメータリカバリ

真の分布が確率モデルに対して正則であり，事後分布が正規分布に近似できる場合，積分のラプラス近似の手法を用いて，自由エネルギーを近似的に求められます．それが，BIC(Bayesian Information Criterion)であり，最尤推定から得られる最大対数尤度を用いた単純な式として，$$\begin{array}{r}\mathrm{BIC}=-\sum _{i=1}^n\log p(x_i|{\hat{w}}_{\mathit{\text{ML}}})+\frac{d}{2}\log n\end{array}$$と定義されます．ただし，$d$はモデルのパラメータ数です．このとき，$$\begin{array}{r}{\mathbb{E}}_{q(X^n)}[\mathrm{BIC}]={\mathbb{E}}_{q(X^n)}[F_n]+O\left(1\right).\end{array}$$が成立します．つまり，サンプルサイズに依存しない1のオーダーの差はあるものの，BICは平均的に自由エネルギーによく近似します．

正則でないモデルにおいても，自由エネルギーの近似値を導出できるのがWBIC(Widely Applicable Bayesian Information Criterion)です．

これは，経験対数損失 $$\begin{array}{r}{L}_n(w)=-\frac{1}{n}\sum _{i=1}^n\log p(x_i|w)\end{array}$$の$n$倍について，逆温度が$\beta =1/\log n$のときの事後分布によって期待値をとったものとして定義されます．つまり，$$\begin{array}{r}\mathrm{WBIC}=\int n{L}_n(w)\left[\frac{\prod _{i=1}^{n}p(x_i|w{)}^{\beta }\phi (w)}{\int \prod _{i=1}^{n}p(x_i|w{)}^{\beta }\phi (w)dw}\right]dw\end{array}$$です．このとき， $$\begin{array}{r}{\mathbb{E}}_{q(X^n)}[\mathrm{WBIC}]={\mathbb{E}}_{q(X^n)}[F_n]+O(\log (\log n)).\end{array}$$が成立します．$\log (\log n)$のオーダーの差は残りますが，WBICは平均的に自由エネルギーに近似します．

そのほか，MCMCサンプルから自由エネルギーを推定する手法としてブリッジ・サンプリング(bridge sampling)があります．

うえで計算した自由エネルギーの厳密値と，WBICやブリッジサンプリングからの計算が一致するかどうかを確認してみましょう．

以下のようなデータを使い，データを生成する未知の分布を確率モデルで推測します．「未知」と言いつつ，ここでは《正解》が分かっているデータをRで作成します．

      n = 500
y <- rnorm(n)
data1 <- list(y=y,n=n)
    

確率モデルと事前分布は
$$\begin{array}{rcl}& y[i]& \sim \mathrm{Normal}(\mu ,1)i=1,2,\dots ,n\\ & \mu & \sim \mathrm{Normal}(0,1)\end{array}$$
です．パラメータの事後分布をMCMCで計算するためのStanコード(fe.stan)は，以下の通りです．

      data{
int n;// no. of respondents real y[n];//
}

parameters {
real mu;
}

model {
target += normal_lpdf(mu|0,1);
target += normal_lpdf(y|mu,1);
}
    

後でブリッジサンプリングを使うために，target記法で書いている点に注意してください．つぎにWBICをMCMCサンプルから計算するためのコード(wbic.stan)です．

      data{
int n;// no. of respondents real y[n];//
}

parameters {
real mu;
}

model {
target += normal_lpdf(mu|0,1);
target += (1/log(n))*normal_lpdf(y|mu,1);
}

generated quantities{
vector[n] log_lik;
for(i in 1:n){
log_lik[i] = normal_lpdf(y[i]|mu,1);
}
    }
    

WBICは逆温度$\beta$が$1/\log n$であるときの事後分布を使った経験対数損失（の$n$倍）の平均値なので，確率モデルを$1/\log n$倍しています．

以上のStanコードを次のようなRコードで実行します．

      library(rstan)
library(bridgesampling)

# define WBIC function 
WBIC.mcmc<-function(x){-mean(rowSums(x))}

############################################################
# compile and sampling
############################################################
model <- stan_model(file='model/fe.stan')
fit <- sampling(model,data=data1,iter = 10500, warmup =500, seed=1234,chains = 5)
bs<-bridge_sampler(fit,method="warp3")
-logml(bs)

model.wbic <- stan_model(file='model/wbic.stan')
fit.wbic <- sampling(model.wbic,data=data1,iter = 5500, warmup =500, seed=1234,chains = 3)
WBIC.mcmc(rstan::extract(fit.wbic)$log_lik)
    

計算した結果です．

      # ブリッジサンプリングから計算した自由エネルギー
# [1] 706.8302, n=500　MCMCsample=40000
# [1] 706.8301, n=500　MCMCsample=50000
# MCMCサンプルから計算したWBIC
# 706.7799, n=500
    

そして最後に解析的に解いた自由エネルギーにデータをあてはめて計算します（ここが今回のノートのハイライトです）．

      a1 <- -sum(log(1/sqrt(2*pi)*exp(-y^2/2)))
nume <- sqrt((2*pi)/(n+1))*exp((sum(y^2))/(2*(n+1)))
a2 <- -log(nume/sqrt(2*pi))
a1+a2
# [1]   707.0134  自由エネルギーの厳密値
    

ほぼ一致しました．ちょびっとズレていますが誤差の範囲でしょう．これで安心してブリッジサンプリングやWBICが使えますね．

文献

• 浜田宏・石田淳・清水裕士, 2019, 『社会科学のためのベイズ統計モデリング』朝倉書店.
• 渡辺澄夫, 2012, 『ベイズ統計の理論と方法』コロナ社.