# 4章　漸近理論と非ベイズアプローチとの関係

## 事後分布の正規近似
以下のノートは、増村ノートにPRML上巻213-214を書き写しています。(齋藤)

まず、$y$がスカラーの場合から考え、ベクトルへの拡張を考える。

一つの最頻値$\hat{\theta}$を持つ事後分布$p(\theta|y)$を想定する．ラプラス近似により、確率密度分布$p(\theta|y)$に対し、正規分布による近似を目指す。まず、 $\log p(\theta|y)$に対して事後分布の最頻値$\hat{\theta}$のまわりでテイラー展開（2次近似）を行うと、

\begin{align}
\log p(\theta|y) &= \log p(\hat{\theta}|y) + (\theta - \hat{\theta})\frac{d}{d\theta}[\log p(\theta|y)]_{\theta = \hat{\theta}} + \frac{1}{2}(\theta - \hat{\theta})^{2}\frac{d^{2}}{d\theta^{2}}[\log p(\theta|y)]_{\theta = \hat{\theta}}　+ \cdots\\
&= \log p(\hat{\theta}|y) + \frac{1}{2}(\theta - \hat{\theta})^{2}\frac{d^{2}}{d\theta^{2}}[\log p(\theta|y)]_{\theta = \hat{\theta}}　+ \cdots \\
&\simeq \log p(\hat{\theta}|y) + \frac{1}{2}(\theta - \hat{\theta})^{2}\frac{d^{2}}{d\theta^{2}}[\log p(\theta|y)]_{\theta = \hat{\theta}}　
\end{align}
となる．今、 仮定より、$\theta = \hat{\theta}$で最頻値を持つ、すなわち極大値をので、 一階微分は0となり、その項は消える。

今、両辺の対数を除くと、
\begin{align}
p(\theta|y)
&\simeq p(\hat{\theta}|y) \exp(\frac{1}{2}(\theta - \hat{\theta})^{2}\frac{d^{2}}{d\theta^{2}}[\log p(\theta|y)]_{\theta = \hat{\theta}}　)\\
&\propto \exp(\frac{1}{2}(\theta - \hat{\theta})^{2}\frac{d^{2}}{d\theta^{2}}[\log p(\theta|y)]_{\theta = \hat{\theta}}　) \\
&= \exp(-\frac{I(\hat{\theta}))}{2}(\theta - \hat{\theta})^2) &\because I(\theta) = - \frac{d^{2}}{d\theta^{2}}\log p(\theta|y))
\end{align}
なお、ガウス積分より、$\exp(-\frac{I(\hat{\theta})}{2}(\theta - \hat{\theta})^2) = \sqrt{\frac{I(\hat{\theta})}{2\pi}}$である。

よって、
\begin{align}
\theta|y \approx N(\hat{\theta}, I(\hat{\theta})^{-1})
\end{align}

同様に、$\mathbf{\theta}$がベクトルの時、上記の式は

\begin{align}
\log p(\mathbf{\theta}|y) &\simeq \log p(\hat{\mathbf{\theta}}|y) + \frac{1}{2}(\mathbf{\theta} - \hat{\mathbf{\theta}})^\top\frac{d^{2}}{d\mathbf{\theta}^{2}}[\log p(\mathbf{\theta}|y)]_{\mathbf{\theta} = \hat{\mathbf{\theta}}}(\mathbf{\theta} - \hat{\mathbf{\theta}})
\end{align}
となる．　同様に、対数を除いて正規分布との対応を目指すと、
\begin{align}
p(\mathbf{\theta}|y)　
&\simeq p(\hat{\mathbf{\theta}}|y) \exp( \frac{1}{2}(\mathbf{\theta} - \hat{\mathbf{\theta}})^\top\frac{d^{2}}{d\mathbf{\theta}^{2}}[\log p(\mathbf{\theta}|y)]_{\mathbf{\theta} = \hat{\mathbf{\theta}}}(\mathbf{\theta} - \hat{\mathbf{\theta}})) \\
&\propto \exp( \frac{1}{2}(\mathbf{\theta} - \hat{\mathbf{\theta}})^\top \mathbf{I(\hat{\theta})}(\mathbf{\theta} - \hat{\mathbf{\theta}})) &(\because \mathbf{I(\theta)} = \frac{d^{2}}{d\mathbf{\theta}^{2}}[\log p(\mathbf{\theta}|y)]_{\mathbf{\theta} = \hat{\mathbf{\theta}}}) &&(4.a)
\end{align}
になる。なお、$\mathbf{I(\hat{\theta})}$はヘッセ行列である。

よって、
\begin{align}
\theta|y \approx N(\hat{\mathbf{\theta}}, [\mathbf{I(\hat{\theta})}]^{-1}) &&(4.2)
\end{align}

**例**：平均、分散ともに未知の正規分布について、$y_{1},y_{2},..., y_{n} \sim N(\mu, \sigma^{2})$, 事前分布として一様分布$(\mu, \log \sigma)$を想定すると、事後分布の対数を取ると、p. 64の(3.2)式より、

\begin{align}
\log p(\mu, \log \sigma | y) = constant - n\log\sigma - \frac{1}{2\sigma^{2}}((n-1)s^{2} + n(\bar{y} - \mu)^{2})
\end{align}
となる．以下、$\mathbf{\theta} = (\hat{\mu}, \log \hat{\sigma}$)、$\mathbf{I(\hat{\theta})} = \mathbf{I}(\hat{\mu}, \log \hat{\sigma})$を導出することにより、(4.a)式を目指す。

この一次導関数は

\begin{align}
\frac{d}{d\mu}\log p(\mu, \log \sigma | y) &= \frac{n(\bar{y} - \mu)}{\sigma^{2}}\\
\frac{d}{d (\log \sigma)}\log p(\mu, \log \sigma | y) &= -n + \frac{(n-1)s^{2} + n(\bar{y} - \mu)^{2}}{\sigma^{2}}
\end{align}

$\because \log \sigma = x$とおくと、
\begin{align}
\frac{d}{d (\log \sigma)}\log p(\mu, \log \sigma | y) &= -n - \frac{1}{2} ((n-1)s^2)+n(\bar{y}-\mu)^2)\frac{d}{dx}\exp(-2x) \\
&= -n - \frac{1}{2} ((n-1)s^2)+n(\bar{y}-\mu)^2)(-2)\exp(-2x) \\
&= -n + ((n-1)s^2)+n(\bar{y}-\mu)^2)\exp(-2x) \\
&= -n + ((n-1)s^2)+n(\bar{y}-\mu)^2)\frac{1}{\sigma^2}
\end{align}

これらが0になるような値を考えると、事後分布におけるそれぞれのパラメーターの最頻値は
\begin{align}
(\hat{\mu}, \log \hat{\sigma}) = (\bar{y}, \log(\sqrt \frac{n-1}{n}s))
\end{align}
となる．

二次導関数は

\begin{align}
\frac{2d}{d\mu^{2}}\log p(\mu, \log \sigma | y) &= -\frac{n}{\sigma^{2}}\\
\frac{d^{2}}{d\mu d(\log \sigma)}\log p(\mu, \log \sigma | y) &= -2n\frac{\bar{y} - \mu}{\sigma^{2}}\\
\frac{d^{2}}{d (\log \sigma)^{2}}\log p(\mu, \log \sigma | y) &= \frac{-2((n-1)s^{2} + n(\bar{y} - \mu)^{2})}{\sigma^{2}}
\end{align}

となる．

最頻値における二次導関数の行列(ヘッセ行列)$\mathbf{I}(\hat{\mu}, \log \hat{\sigma})$は、$\hat{\mu} = \bar{y}, \hat{\sigma} = \sqrt \frac{n-1}{n}s$を代入することで以下のようになる．

\begin{pmatrix}
\frac{-n}{\hat{\sigma}^{2}} & 0 \\
0 & -2n \\
\end{pmatrix}

よって、正規分布の近似は

\begin{align}
\log p(\mu, \log \sigma | y) \approx N(\begin{pmatrix}
\mu \\
\log \sigma
\end{pmatrix}|\begin{pmatrix}
\bar{y} \\
\log \hat{\sigma}
\end{pmatrix}, \begin{pmatrix}
\frac{\hat{\sigma}^{2}}{n} & 0 \\
0 & \frac{1}{2n} \\
\end{pmatrix})
\end{align}

- $d$ dimentional の正規分布において、$exp(\chi^{2}_{d}の95th\ percentile/2)$が事後分布の95%点に対応する．
- 正規近似がいつ適当なのかに関する一般的なガイドラインを示すのは難しい．　
- $\theta$の変形によって収束は向上する．
- $\hat{\theta}, I(\hat{\theta})$は十分統計量．
- nが有限の時、full joint dist.よりもconditional dist.やmarginal dist.の方が正確なことが多い．
- $\theta$が低次元のとき、正規近似は頻繁に受け入れられる．高次元の時は、個々の$\theta$に関しては正規近似が使えることが多く、また、正規近似できそうな$\theta$と出来なさそうな$\theta$に分けて、conditional dist.をとることは有用．

## **4.2** 大数の法則
asymptotic normality: データが多くなればなるほど、事後分布は多項正規分布に近づく．

**Consistency（一致性）**: $n \rightarrow \infty$になるにつれ推定量が真の値に近づく性質．
$f(y)$ を真のデータ分布、このデータがあるパラメトリック族（parametric family）$p(y|\theta)$と事前分布$p(\theta)$でモデル化されているとき、

- $f(y) = p(y|\theta_{0})$がある$\theta_{0}$で成立しているとき（＝データ分布はそのパラメトリック族に含まれているとき）、Consistencyがあり、事後分布は$\theta_{0}$に近づく．
- $f(y) = p(y|\theta_{0})$がある$\theta_{0}$で成立してい**ない**とき、 $\theta_{0}$を$f(\theta)$と$p(y|\theta_{0})$の距離が最小になる$\theta$（Kullback-Leibler divergence（カルバック・ライブラー情報量）を最小にする$\theta$）とおくと、事後分布は平均$\theta_{0}$分散$(nJ(\theta_{0}))^{-1}$の正規分布になる．　（ここで$J()$はフィッシャー情報量）

また、$\hat{\theta}$を事後分布$p(\theta|y)$の最頻値とすると、$\hat{\theta} \rightarrow \theta_{0}$、as $n \rightarrow \infty$.

ここで最初にテイラー展開した式の二次の項の係数は以下のように書き換えることができる．
\begin{align}
&p(\theta|y) \propto p(\theta)p(\mathbf{y}|\theta) \\
&\Leftrightarrow \log p(\theta|\mathbf{y}) = \log p(\theta) +\sum_i \log p(y_i|\theta) + const. \\
&\Leftrightarrow \frac{d^{2}}{d\theta^{2}} \log p(\theta|\mathbf{y}) = \frac{d^{2}}{d\theta^{2}} \log p(\theta) + \frac{d^{2}}{d\theta^{2}} \sum_i \log p(y_i|\theta) \\
&\Rightarrow [\frac{d^{2}}{d\theta^{2}}\log p(\theta|\mathbf{y})]_{\theta = \hat{\theta}} = [\frac{d^{2}}{d\theta^{2}}\log p(\theta)]_{\theta = \hat{\theta}} + \sum^{n}_{1}[\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}}
\end{align}

これは定数＋n個の尤度の項の合計として捉えられる． 尤度の個々の平均は、$f(y) = p(y|\theta_{0})$がある$\theta_{0}$で成立しているときは$-J(\theta_{0})$、そうでない時は$E_{f}(\frac{d^{2}}{d\theta^{2}}\log p(y|\theta))$を$\theta=\theta_{0}$で評価した値となる．　これらの値は負であるため、$[\frac{d^{2}}{d\theta^{2}}\log p(\theta|y)]_{\theta = \hat{\theta}}$はnが増えるに従って増加する（？）．

齋藤メモ p. 589のConsiderから始まる段落の下から5行目''is thus negative''は''positive''の間違いでは?KL divergenceの最小値の話をしており、最小値の2階微分は正でないとおかしい、そう考えると、$\sum^{n}_{1}[\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}} = nE_f([\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}}) >0$であり、'increases with order n'に合致する。

\begin{align}
\lim_{n \rightarrow \infty, \hat{\theta} \rightarrow \theta_0} [\frac{d^{2}}{d\theta^{2}}\log p(\theta|\mathbf{y})]_{\theta = \hat{\theta}}
&= \lim_{n \rightarrow \infty, \hat{\theta} \rightarrow \theta_0} ([\frac{d^{2}}{d\theta^{2}}\log p(\theta)]_{\theta = \hat{\theta}} + \sum^{n}_{1}[\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}} )\\
&= nE([\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \theta_0})) \\
&= -nJ(\theta_0) &(\because (2.20)
\end{align}

中村メモ
ここでは教科書p.589の記述、
...the mean is $E_f([\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}})$, which is the negative second derivative of the Kullback-Leibler divergence, $KL(\theta_0)$, and is thus negative, because $\theta_0$ is defined as the point at which $KL(\theta)$ is minimized.
について考えたことを記す。

まず、KLダイバージェンスは以下のように定義される。
\begin{align}
KL(\theta) &= E[\log\frac{f(y_i)}{p(y_i|\theta))}]\\
&= E[\log f(y_i)] - E[\log p(y_i|\theta)]
\end{align}
ここで、本文中の(negative)second termとは$-E[\log p(y_i|\theta)]$のこと。
まず、KLダイバージェンスを$\theta$で２回微分すると第1項が消え、第２項の2回微分、$-E_f([\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}})$になることがまずわかる。
今、真のモデルパラメータ$\theta_0$について考えると、$KL(\theta_0)$の第1項は定数と同様に考えられる($\theta$を含んでいない)ので、KLダイバージェンスを最小化させることはnegative second term、つまり$-E[\log p(y_i|\theta_0)]$を最小化させることに相当する。従って、negative second termは負の値となる。これはすなわち$E[\log p(y_i|\theta_0)]$の最大化である。従って、テイラー展開における$E_f([\frac{d^{2}}{d\theta^{2}}\log p(y_{i}|\theta)]_{\theta = \hat{\theta}})$は上の齊藤さんの議論からもnのオーダーで増えていくことがわかる。

まとめ．1　$\hat{\theta} \rightarrow \theta_{0}$、as $n \rightarrow \infty$.
　　　　2　テイラー展開の第二項の係数の負の部分は$-nJ(\theta_{0})$に近づく．
　　　　3 nが増えるに従い、尤度関数が事前分布より優位になるため、正規近似をする際にはモデ
　　　　　ルの尤度があれば十分になる．

## 4.3 漸近理論への反例
- モデルがunderidentifiedなとき
尤度　$p(y|\theta)$が　$\theta$の範囲と同じとき、$p(y|\theta)$はフラットな尤度と呼ばれ、収束する先の一つの$\theta_{0}$が存在しない．　（データがパラメーターに関する情報を持ってないため、nの大きさに関わらず事後分布のパラメーターは事前分布と同じになってしまう）　

- パラメーターの数がnとともに増加するとき
　例えば、$y_{i} \sim N(\theta_{i}, \sigma^{2})$の時、ユニット毎に取れるデータがｎが増えるごとに増えない限り、パラメーター$\theta_{i}$は一致性をもって推定することはできない。

- aliasing (エイリアシング)：underidentifiedの特殊ケース
例）　$p(y_{i}|\mu_{1},\mu_{2}, \sigma^{2}_{1},\sigma^{2}_{2}, \lambda) = \lambda\frac{1}{\sqrt{2\pi}\sigma_{1}}e^{-\frac{1}{2\sigma^{2}_{1}}(y_{i}-\mu_{1})^{2}} + (1-\lambda)\frac{1}{\sqrt{2\pi}\sigma_{2}}e^{-\frac{1}{2\sigma^{2}_{2}}(y_{i}-\mu_{2})^{2}}$ というモデルを想定する。
ここで、$(\mu_{1},\mu_{2})と( \sigma^{2}_{1},\sigma^{2}_{2})$や、$\lambda$と$(1-\lambda)$をいれかえても変わらない。そのため、一点に収束しない。
この問題はパラメーター空間を制限し、重複が起きないようにすれば解決できる。

- 非有界の尤度
尤度関数が有界でない場合、事後分布がパラメーター空間にない可能性があり、一致性も正規近似も成立しない。e.g.）ｎが増えるに従って、尤度の値が無限に近づくときなど
このような状況は実際にはほとんど見られない。もし起きてもいくつか解決策がある（ｐ９０参照）

- 不適切な事前分布
事後分布を計算する際に、積分が無限になってしまうような事前分布は不適切。（？）
不適切な事前分布を設定すると、正規近似は成立しない。
例）ｎ回成功、０回失敗のバイノミアル分布に対して事前分布$Beta(0,0)$を適用する。
解決法としては、そもそも不適切な事前分布を用いないか、用いたとしても有限な積分になる事前分布を用いる。

- 収束点を除外している事前分布
もし、$p(\theta_{0}) = 0$なら漸近理論の結果は成立しない。
少しでも収束点になる可能性がある$\theta$全てに正の確率をあたえるような事前分布を用いる。

- パラメータ空間の端に収束する場合
もし、$p(\theta_{0})$がパラメーター空間上の端にある場合、極限においても正規近似は適当でない
少しでも収束点になる可能性がある$\theta$全てに正の確率をあたえるような事前分布を用いる。

- 分布の端
正規近似した分布の両端では、正規近似は正確ではない。

## 4.4 ベイズ推論の頻度論的評価
- Large sample correspondence
事後分布$p(\theta|y)$の正規近似は．
\begin{align}
p(\theta|y) \approx N(\hat{\theta}, [I(\hat{\theta})]^{-1})
\end{align}

, where $I(\theta) = - \frac{d^{2}}{d\theta^{2}}\log p(\theta|y)$.
であわらせられた。

これを変形して

\begin{align}
[I(\hat{\theta})]^{1/2}(\theta - \hat{\theta})|y \sim N(0, I)
\end{align}
, where $\hat{\theta}$ が事後分布の最頻値で、$[I(\hat{\theta})]^{1/2}$が行列$I(\hat{\theta})$の平方根

（参考：　$X \sim N(\mu, \sigma^{2})$ならば$\frac{X-\mu}{\sigma} \sim N(0, 1)$）

もし、$f(y) \equiv p(y|\theta)$がある$\theta$について成立するならば、固定された$\theta$に対して$n \rightarrow \infty$において
\begin{align}
[I(\hat{\theta})]^{1/2}(\theta - \hat{\theta})|\theta \sim N(0, I)
\end{align}
が成立する。（？）

これは、$(\theta-\hat{\theta})$のどんな関数に対しても、２つ目の式から生じた事後分布は近似的に３つ目の式からの繰り返しサンプリングした結果と同じであることを示している。よって、事後分布の９５％信用区間は繰り返しサンプリングの９５％が真の値を含んでいる。

- 点推定、一致性、有効性
ｎが大きいときは「点推定」の結果はある程度信用できるが、ｎが小さなときは点ではなく事後分布そのものを見るべき。
$f(y) \equiv p(y|\theta)$が成立していたら、事後分布の最頻値$\hat{\theta}$だけでなく、平均や中央値なども、一致性（もしくは、近似的に不偏、asymptotic unbiassness。一致性と似た概念）と有効性（efficiency）をもつ

（参考：有効性とは、他のパラメーターの推定量との比較において、その推定量の分散がもっとも小さいことを指す）

- 信用区間
下手に意訳してまちがえるより教科書の説明をそのまま載せたほうが良さそうなので、コピペします。
$100(1-\alpha)\%$ confidence regionとは、
If one chooses $\alpha$ to be small enough (for example, 0.05 or 0.01), then since confidence regions cover the truth in at least (1 − $\alpha$) of their applications, one should be confident in each application that the truth is within the region and therefore act as if it is.

## 4.5 他の統計手法のベイズ的解釈

- ある確率モデルから生じる大きなサンプルを扱う問題では、ベイズは他の統計手法と似ている。
- ｎが小さいときでも、他の手法はベイズの近似として捉えることが可能である。その際、impricitな事前分布を考えることは有用

- 最尤推定
極限において、最尤推定量$\hat{\theta}$は十分統計量。→事後分布の最頻値、平均、中央値も十分統計量。ｎが増えるごとに事前分布は関係なくなっていくので、都合の良い無情報事前分布を用いることが正当化される。
ただ、ｎが限られてくる状況では、$\hat{\theta}$は十分統計量ではなく、そのため最尤推定は不十分。パラメーターの数が多いときでは無情報事前分布は正当化しづらく、一致性は有用ではない。そのようなときはCh.5での階層モデルを用いる。

- 不偏推定量
ベイジアンの視点からは不偏推定量はｎが大きいときは妥当性があるが、ｎが少ないときはそうではない。
不偏推定量の問題は複数のパラメーターを一度に求めるのが不可能なこと、そして将来の観察可能な値を予測問題におけるパラメータとして扱っていること。また、未知の値を予測と推定のどちらかに区分しなければいけないが、両者の区別は曖昧。

- 信頼区間
信用区間と似ているときもあるが、信頼区間が０や無限大になったりすることや、現実的な範囲の９５％をカバーしている信頼区間があったとき、繰り返しサンプリングしたものの９５％に真の値を含んでいることは当たり前だったりと、いくつか不可解な点がある。

- 仮説検定
対抗仮設と帰無仮説に分ける二分法は解釈が難しいし、本来は帰無仮説のパラメーターが０かどうかではなく、事後分布に興味がある。仮説検定はモデルの適合度を測るのに有用かもしれない。

- マルチレベル　モデリング
複数の比較を含む問題を扱う際には、階層モデルがよい。複数の比較を含む問題の懸念点はベイズで自動的に対処できる

- ノンパラメトリック手法、Permutation test（並び替え検定）、ジャックナイフ、ブートストラップ
順位和検定での仮説検定などベイズの枠組みからは解釈しにくいものが多い。ジャックナイフ、ブートストラップなども具体的な確率モデルがないと解釈しにくい。ベイズの観点からは同時分布（joint prob.）を作り、データに合致するか見たほうが良い。
Permutation testなどいくつかのノンパラメトリックな手法は無情報事前分布を用いたベイズ統計と似た結果になる。しかし、複数のパラメーターや説明変数、有情報事前分布を入れるとうまく行かない。
例えば、ウィルコクソンの順位和検定はベイズの近似として捉えることができ、そのような見方をすることで独立変数や欠損データを扱えるし、model-based approachの一般性も示している。