尤度とAIC

KL情報量

データの当てはまりのよさを数値で表すことを考える． これを実現するために Kullbuck と Leibler が提案した量を Kullbuck-Leibler 情報量または KL 情報量と呼ぶ． 離散型確率変数の場合， ${\boldsymbol p}=\{p_1, p_2, \cdots, p_n\}$ を 真の確率分布， ${\boldsymbol q}=\{q_1, q_2, \cdots, q_n\}$ をモデルの確率分布とおくと，

$$I({\boldsymbol p}, {\boldsymbol q})\equiv \sum_{i=1}^{n} p_i \log p_i/q_i$$ (51)

ここで，$n$ は確率変数の数であり， $\log$ は自然対数を表す． 連続型確率変数の場合は，真の分布を $p(x)$ モデルの分布を $q(x)$ とすると，

$$I(p, q)\equiv \int_{-\infty}^{\infty} p(x) \log p(x)/q(x) dx$$ (52)

である． KL情報量は以下のような性質を持っている．

1. KL情報量は常に0以上である．
2. モデルの分布が真の分布と一致するときにKL情報量はゼロとなる．

この性質はKL情報量をモデルの当てはまりのよさを与えるひとつのものさしとす ることを支持している．また，KL情報量が（負の）エントロピーに符号を逆にし たものであることも，ものさしとすることを支持するものである． 上の性質から明らかなように，このものさしはあてはまりが良いほど小さい値を 持つ． 以下では，このものさしを用いてモデルの当てはまりのよさを比較する方法につ いて考えてゆく．

対数尤度

KL情報量は，離散的な場合，

$$I({\boldsymbol p}, {\boldsymbol q}) \equiv \sum_{i=1}^{n} p_i \log p_i/q_i$$

$$= \sum_{i=1}^{n} p_i \log p_i-\sum_{i=1}^n q_i \log q_i$$ (53)

$$= E[\log p]-E[\log q]$$ (54)

連続的な場合，

$$I(p, q) \equiv \int_{0}^{\infty} p(x) \log p(x)/q(x) dx$$

$$= \int_{-\infty}^{\infty} p(x)\log p(x) -\int_{-\infty}^\infty p(x)\log q(x) dx$$

$$= E[\log p]-E[\log q]$$ (55)

と変形できる．どちらの式も，第一項は真の分布から得られる定数であり， モデルの当てはまりのよさを決めるのは第二項であることがわかる． この第二項の符号を逆転した式 $E[\log q]$ の近似値をデータから得ることを 考える． 観測データを ${\boldsymbol x} =\{x_1, x_2, \cdots, x_m\}$ とする． 離散型確率変数の場合で $x_i$ が出現する確率を $q(x_i)$ ， 連続型確率変数の場合で $x_i$ についての確率密度関数をやはり $q(x_i)$ と 書くと，上記の $E[\log q]$ の近似値は，

$$\frac{1}{m} \sum_{i=1}^n \log q(x_i)$$ (56)

と表せそうである． これを平均対数尤度と呼ぶ．これをデータ数 $m$ 倍した

$$\ell \equiv \sum_{i=1}^n \log q(x_i)$$ (57)

を対数尤度（log likelihood）と呼ぶ．また，対数をとる前の

$$\prod_{i=1}^n q(x_i)$$ (58)

を尤度（likelihood）と呼ぶ．Likelihood（もっともらしさの度合い）というの は，この式が，考えている確率モデルにおけるそれぞれの観測データの積事象が 観測される確率になっていることからも理解できる． 対数尤度や尤度は確率モデルのよさを測るものさしとなる． 対数尤度とKL情報量を比較してみると符号が逆転しているので，KL情報量とは逆に， モデルの対数尤度が大きいほど，当てはまりがよくなることに注意する．

最尤法

定義

考えているモデルの尤度が最大になるように パラメータの値を決めることによって，モデルが真の分布に 最も近くなるように調整することを最尤法という． また，最尤法によって得られたパラメータの値を最尤推定量と呼び， 変数名の上にハット（＾）をつけて表す．

正規分布

$n$ 個のデータ $x_1, x_2, \dots, x_n$ が正規分布する，すなわち $x$ の確率密度関数が，

$$f(x\vert\mu , \sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left\{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right\}$$ (59)

で表される仮定できるとする． このときにデータから得られる最尤推定量 $\hat{\mu}$ および $\hat{\sigma}^2$ についての式を求めることを考える． このとき，対数尤度は，

$$\ell(\mu, \sigma^2) = \sum_{i=1}^n \log f(x_i\vert\mu, \sigma^2)$$

$$= \sum_{i=1}^n \left\{\log \sqrt{\frac{1}{2\pi \sigma^2}}+ \frac{1}{2}\left(\frac{x_i-\mu}{\sigma}\right)^2\right\}$$

$$= -\frac{n}{2} \log (2\pi \sigma^2)-\frac{1}{2\sigma^2} \sum_{i=1}^n (x_i-\mu)^2$$ (60)

となる．この量が最大となるのは，

$$\frac{\partial \ell}{\partial \mu}=0, \qquad \frac{\partial \ell}{\partial \sigma^2}=0$$ (61)

のとき．この方程式の解 $\hat{\mu}$ と $\hat{\sigma}^2$ は，

$$\hat{\mu} = \frac{1}{n}\sum_{i=1}^n x_i$$ (62)

$$\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i-\hat{\mu})^2$$ (63)

となる．すなわち，あるデータセットについて，それが正規分布すると 仮定して最尤法によって推定されたパラメータの値は，平均の推定値が平均値， 分散の推定値が標本分散になることを示している． 後のために必要になるので，このときの対数尤度 $\ell(\hat \mu, \hat \sigma^2)$ を求めておくと，

$$\ell(\hat \mu, \hat \sigma^2)=-\frac{n}{2}\log (2\pi\hat\sigma^2)-\frac{n}{2}$$ (64)

となることがわかる．

多項式回帰モデル

$n$ データのペア $(x_1, y_1), (x_2, y_2), \cdots (x_n, y_n)$ を 考える． これらのデータが，

$$y=a_0 +a_1 x + a_2 x^2 + \cdots + a_{m-1} x^{m-1} + a_m x^m$$ (65)

に従うと仮定し，このときの AIC を求めることを考える． そのために，実際に得られたデータは平均ゼロ・分散 $\sigma^2$ の 誤差を伴うとしたした確率モデルを考える． このようなモデルを$m$次多項式回帰モデルという． このモデルの対数尤度は，誤差が正規分布に従うという 確率モデルであることから，

$$\log q_i = \log \left[ \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left\{ \frac{\left( y_i- a_0 -a_1 x_i - \cdots - a_m x_i^m \right)^2} {2\sigma^2} \right\} \right]$$

$$= -\frac{1}{2} \log 2\pi -\frac{1}{2}\log \sigma^2 -\frac{1}{2\sigma^2} \left( y_i- a_0 -a_1 x_i - \cdots - a_m x_i^m \right)^2$$ (66)

と表されることにより，

$$\ell(a_0, a_1, \cdots , a_m, \sigma^2) = \sum_{i=1}^n \log q_i$$ (67)

$$= -\frac{n}{2} \log 2\pi -\frac{n}{2}\log \sigma^2$$

$$-\frac{1}{2\sigma^2}\sum_{i=1}^n \left\{ y_i -\left( a_0 +a_1 x_i + \cdots + a_m x_i^m \right)\right\}^2$$ (68)

このとき，最尤推定量（$\hat{a}_0$, $\hat a_1$, $\cdots$, $\hat a_m$, $\hat\sigma^2$ ）は

$$\frac{\partial l}{\partial a_0}= \frac{\partial l}{\partial a_1}... ...dots =\frac{\partial l}{\partial a_m}= \frac{\partial l}{\partial \sigma^2}=0$$ (69)

の解である．この方程式は最初の $(m+1)$ 個について以下のように書き換えられる．

$$\begin{bmatrix} n & \sum x_i & \sum x_i^2 & \cdots & \sum x_i^m ... ...in{bmatrix} \sum y_i\\ \sum x_i y_i \\ \vdots \\ \sum x_i^m y_i \end{bmatrix}$$ (70)

であり，この解は最小自乗法によって得られるものと同じである． また，最後の方程式は，

$$-\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2} \sum_{i=1}^{n}(y_i-\hat a_0-\hat a_1 x_i -\cdots -\hat a _m x_i^m)^2=0$$ (71)

の解であって，

$$\hat{\sigma}^2=\frac{1}{n} \sum_{i=1}^{n}(y_i-\hat a_0-\hat a_1 x_i -\cdots -\hat a _m x_i^m)^2$$ (72)

である． 後のために必要になるので，このときの対数尤度を求めておくと，

$$\ell(\hat{a}_0, \hat{a}_1, \cdots ,\hat{a}_m, \hat{\sigma}^2) = ... ...^2 -\frac{n}{2}= -\frac{n}{2} \left(\log 2\pi + 1 +\log \hat{\sigma}^2\right)$$ (73)

AIC

AICとは

$n$ 次回帰モデルの最尤法の結果について改めて考えてみると， ある問題点に気が付く． それは次数が大きいほど残差平方和は小さくなることである． つまり，残差平方和が小さいことが正義だとすると，次数を 大きくとることが正しいモデル選択の方法になってしまう． 極論を言うと $n$ 組のデータを用いて $n$ 次回帰モデルを 作ることが一番よいモデルであるということになる． しかし，この考え方はどうみてもおかしい． 赤池の情報量規準(Akaike's Information Criterion: AIC)は， このような問題点を回避するためのものさしとなるものである．

AICの導出

AICの導出については参考図書を見てほしい．ここでは， 何をもとにしてどんな仮定を行うことで導出されるのかを 見るに留めておく． 基本となるのは KL 情報量である．これを以下の条件（仮定）の もとで評価してみる．

• 最尤モデルの平均対数尤度を真のパラメータで Taylor 展開したときに２ 次の項まで考慮すれば十分である．
• 最尤推定量はデータの個数が無限大になるときに，真のパラメータとなる．

これらを仮定すると，パラメータ数 $K$ の期待対数尤度 （平均対数尤度のデータ ${\boldsymbol x}$ に関する期待値） $l_n^*(K)=$ つまり，ＫＬ情報量

$$I(p,q) = E[\log (p)]-E[\log (q)]$$ (74)

の第２項が，

$$\frac{(最大対数尤度)-(パラメータ数)}{n}$$ (75)

で近似されることがわかる（$n$ はデータ数)． そこで，（歴史的な経緯から）分子に (-2) を 乗じたものを ${\rm AIC}(K)$ と呼び，この数値によって モデルのあてはまりの良さを評価する． すなわち，モデルのパラメータ数を $K$ としたとき，

$${\rm AIC(K)}=-2\ell(\cdots)+2K$$ (76)

が小さいモデルほど良いモデルであるとして評価を行う．

AICの使用例

硬貨の真偽

$n$ 回の試行で $r$ 回が表であった．このとき，この硬貨の真偽を判断してみ る．ふたつのモデルが考えられる．ひとつは硬貨が本物であると考えるモデルで ある．これをモデルＡと呼ぶことにする．モデルＡでは表がでる確率は 0.5 と みなされる．もうひとつは硬貨が偽物であると考えるモデルである． これをモデルＢと呼ぶことにする．モデルＢでは表がでる確率 $p$ が不明なので， 最尤法を用いて $p$ の値を求める． 表がでる確率が $p$ であるとき，その対数尤度 $l(p)$ は，

$$l(p)=\sum_{i=1}^{n} \log p_i =\sum_{i=1}^r \log p + \sum_{i=1}^{n-r} \log(1-p) =r\log p+(n-r)\log (1-p)$$ (77)

$l(p)$ が最大になるように最尤推定量 $\hat p$ を設定するには，

$$\frac{\partial l}{\partial p}= \frac{r}{p}-\frac{n-r}{1-p}=0$$ (78)

を解いて，

$$\hat p=\frac{r}{n}$$ (79)

このとき，最大対数尤度は，

$$l(\hat p)=r\log \frac{r}{n}+(n-r)\log \frac{n-r}{n}$$ (80)

である． この結果，モデルＡの AIC は，

$$\textrm{AIC}(0)=-2 \{r \log 0.5 + (n-r) \log 0.5\}+2\times 0 = -2 n \log 0.5$$ (81)

モデルＢの AIC は，

$$\textrm{AIC}(1) = -2 \left\{ r\log \frac{r}{n}+ (n-r)\log \frac{n-r}{n}\right\} +2 \times 1$$

$$= 2\left\{1-r\log \frac{r}{n}-(n-r)\log \frac{n-r}{n} \right\}$$ (82)

考え方は逆になってしまうが， $\hat p$ を固定して，各試行回数に おける AIC を計算した結果を下の表に示した． 100回コインを投げて60回表になったら偽物と思った方がよさそうである．

0.55			0.6
n	AIC(0)	AIC(1)	n	AIC(0)	AIC(1)
50	69.3	70.8	50	69.3	69.3
100	138.6	139.6	100	138.6	136.6
200	277.3	277.3	200	277.3	271.2
300	415.9	414.9	300	415.9	405.8
400	554.5	552.5	400	554.5	540.4

正規分布モデルのAIC

正規分布に関する確率モデルとしては以下の４つが考えられる．

• データはもともと与えられた $\mu$ と $\sigma^2$ に従う．
• 平均は最尤推定量 $\hat \mu$ が正しい．
• 分散は最尤推定量 $\hat\sigma^2$ が正しい．
• 平均，分散ともに最尤推定量 $\hat \mu$ ， $\hat\sigma^2$ が正し い．

これら４つのモデルをそれぞれ，モデルＡ，モデルＢ，モデルＣ，そしてモデル Ｄとする．これらの AIC を $\textrm{AIC}(0)$ ， $\textrm{AIC}(1)$ ， $\textrm{AIC'}(1)$ ， $\textrm{AIC}(2)$ とし， 3.3.2節の結果を用いて求めることにする． モデルＡの AIC は， 対数尤度の式(61)を用いて，

$$\textrm{AIC}(0)=-2\times l(\mu, \sigma^2)+2 \times 0 = n\log (2\pi \sigma^2)+\frac{1}{\sigma^2} \sum_{i=1}^n (x_i-\mu)^2$$ (83)

となる．モデルＢ，モデルＣの結果は一部に最尤推定量を使うので，それぞれ，

$$\textrm{AIC}(1) = -2\times l(\hat \mu, \sigma^2)+2 \times$$

$$= n\log (2\pi \sigma^2)+\frac{1}{\sigma^2} \sum_{i=1}^n (x_i-\hat \mu)^2 +2$$

$$= n\log (2\pi \sigma^2)+n\frac{\hat \sigma ^2}{\sigma^2}+2$$ (84)

$$\textrm{AIC'}(1)=-2\times l(\mu, \hat \sigma^2)+2 \times 1 = \log (2\pi \hat \sigma^2)+\frac{1}{\hat \sigma^2} \sum_{i=1}^n (x_i-\mu)^2$$ (85)

となる．モデルＤはすべて最尤推定量で構成されるので，最大対数尤度の式であ る式(65) を用いて，

$$\textrm{AIC}(2)=-2\times l(\hat \mu, \hat \sigma^2)+2 \times 2 = n\log (2\pi\hat\sigma^2)+n +4$$ (86)

となる．例えば機械の故障の判定などは，これらの４つのモデルを比較すること によって結論を得ることができる．

多項式回帰モデルのAIC

3.3.3節の結果より， $n$ データのペア $(x_1, y_1), (x_2, y_2), \cdots (x_n, y_n)$ を $m$ 次回帰モデル

$$y=a_0 +a_1 x + a_2 x^2 + \cdots + a_{m-1} x^{m-1} + a_m x^m$$ (87)

に当てはめたときの AIC は，

$$\textrm{AIC}(m) = -2 l(\hat{a}_0, \hat{a}_1, \cdots ,\hat{a}_m, \hat{\sigma}^2) +2 \times (m+2)$$

$$= n\left(\log 2\pi + 1 +\log \hat{\sigma}^2\right) +2(m+2)$$ (88)

となる． 例えば，

$i$	1	2	3	4	5	6	7	8	9	10	11
$x_i$	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
$y_i$	0.012	0.121	-0.0194	-0.0183	-0.032	-0.037	0.196	0.077	0.343	0.448	0.434

というデータの組について，多項式回帰モデルを当てはめる． 下の表に各次数における最尤推定量とそのときの最大対数尤度およびＡＩＣをま とめた．パラメータの数を多くすれば多くするほど分散の最尤推定量 （残差平方和/データ数）が小さくなるので，対数尤度が大きくなり モデルを比較するのが困難であること．そして，ＡＩＣを用いると その比較ができ，今回の場合，２次モデルが最良のモデルであると 予想できることがわかる． 実際にグラフに書いてみると２次と３次では曲線の形にほとんど違いがなく， パラメータの数を多くしてまで３次曲線で近似する必要がみられないことがわか る．

モデル	$\hat a_0$	$\hat a_1$	$\hat a_2$	$\hat a_3$	最大対数尤度	AIC
０次	0.138573				3.23417	-2.46835
１次	-0.0852364	0.447618			8.50008	-11.00017
２次	0.0602566	-0.522335	0.969953		13.37498	-18.74997
３次	0.078849	-0.817746	1.74463	-0.516453	13.53748	-17.07496

ARモデルのAIC

AR モデルとは，自己回帰モデルとも呼ばれ，等時間間隔で 計測された時系列に対して，現在の自分を過去の自分の線形和と して表すものである． これは過去に得られたデータ $x_1, x_2, \dots, x_n$ をもとにして， 現在の自分の値を

$$x_t = a_1 x_{t-1} + a_2 x_{t-2} \cdots +a_m x_{t-m}$$ (89)

$$= \sum_{i=1}^{m} a_m x_{t-m}$$ (90)

と近似して表そうというものである． このモデルには平均ゼロ，分散 $\sigma^2$ である 誤差が伴うという確率モデルについて検討する． 知りたいのは最適な $m$ の値とそのときの最尤推定量$\hat a_i$ $(i=1,2,\dots,m)$ の値である． 最初にしなければならないのはデータの番号付けのやり直しである． もとのデータ $x_1, x_2, \dots, x_n$ を $x_{1-m'}, x_{2-m'}, \dots,x_0,x_1,x_2,\dots , x_{n-m'}$ と 置きなおす．このときの $m'$ とは，比較検討する $m'$ の 値の最大値である． $x_{1-m'}, x_{2-m'}, \dots,x_0$ というデータは 観測データの個数には入れない．すなわち， $x_1, x_2, \dots, x_{n-m'}$ を観測データとする．これは $x_1$ の推定値を計算するのに， それよりも過去のデータが必要なためである．以下ではデータの個数として， $n'=n-m'$ を用いて議論する． 誤差が正規分布すると仮定しているので，このモデルの対数尤度は，

$$l(a_1, a_2, \dots, a_m, \sigma^2)= -\frac{n}{2}\log (2\pi \sigma^... ...rac{1}{2\sigma^2} \sum_{t=1}^{n'} \left(x_t-\sum_{i=1}^m a_i x_{t-i}\right)^2$$ (91)

この式はデータとして $x_{1-m}$ までを参照する．前述のようにデータの番号 付けを変更する意味がここでわかる筈である． 最尤推定量を求める方程式は，

$$\frac{\partial l}{\partial a_1}= \frac{\partial l}{\partial a_2}... ...dots =\frac{\partial l}{\partial a_m}= \frac{\partial l}{\partial \sigma^2}=0$$ (92)

となり，最後の方程式以外は，

$$\begin{bmatrix} C(1,1) & C(1,2) & C(1,3) & \cdots & C(1,m) \\ ... ...x} = \begin{bmatrix} C(1,0) \\ C(2,0) \\ \vdots \\ C(m,0) \\ \end{bmatrix}$$ (93)

とまとめられる．ここで $C(i,j)$ とは，

$$C(i,j)=\sum_{t=1}^{n'}x_{t-i}x_{t-j}$$ (94)

である．また最後の方程式を解くと，最尤推定量 $\hat\sigma^2$ を 求めることができ，その結果は，

$$\hat \sigma^2 = \frac{1}{n'}\sum_{t=1}^{n'} \left(x_t-\sum_{i=1... ..._{t-i}\right)^2 =\frac{1}{n'}\left(C(0,0)-\sum_{i=1}^m \hat a_{i}C(i,0)\right)$$ (95)

となる．