AICとBICの導出を理解したい(前編)

#はじめに

- 渡辺澄夫 [データ解析(2017)](http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/da2017.html) 「(13)モデルの評価」「(14)モデルの評価」の読解メモ。
- AICとBICの導出について理解できる範囲でまとめた。

#AICの導出

- 真の分布を$q(x)$、予測モデルを$p(x|w)$とする。$w$は予測モデルのパラメータ。

```math
\begin{align}
G(w)&=-\int q(x)\log p(x|w)dx\quad 汎化損失\\
T(w)&=-\frac{1}{n}\sum_{i=1}^n \log p(x_i|w)\quad 学習損失
\end{align}
```

- 学習損失$T(w)$を最小化するパラメータを$\hat{w}$とする。$\hat{w}$は最尤推定量と呼ばれる。
- 本来は$G(w)$を最小化するパラメータを求めたいが、真の分布$q(x)$が不明なので計算できない。
- そこで$T(\hat{w})$から$G(\hat{w})$を推定する方法を考える。
- 真の分布$q(x)=p(x|w_0)$と表現できる場合を考える。
- ある条件の下で最尤推定量の漸近正規性が成立する。$I(w)$はフィッシャー情報行列。

```math
\begin{align}
\sqrt{n}(\hat{w}-w_0)&\sim N(0,I(w_0)^{-1})\\
I_{jk}(w)&=\int \left[\frac{\partial}{\partial w_j}\log p(x|w)\right]
\left[\frac{\partial}{\partial w_k}\log p(x|w)\right]p(x|w)\mathrm{d}x\\
&=-\int\left[\frac{\partial^2}{\partial w_j\partial w_k}\log p(x|w)\right]p(x|w)\mathrm{d}x
\end{align}
```

- 平均的に以下の関係式が成立する。$d$はパラメータ数。

```math
\begin{align}
G(w_0)-T(w_0)&\simeq 0\quad (q(x)=p(x|w_0))\\
G(\hat{w})-G(w_0)&\simeq\frac{d}{2n}\quad(Gはw_0で最小値)\\
T(w_0)-T(\hat{w})&\simeq\frac{d}{2n}\quad(T は\hat{w}で最小値)
\end{align}
```

- まとめると汎化損失$G(\hat{w})$は学習損失$T(\hat{w})$よりも平均的に$\frac{d}{n}$だけ大きい。
- $2nG(\hat{w})\simeq 2nT(\hat{w})+2d\Rightarrow \mathrm{AIC}$

#正規分布の尤度関数

- 最尤推定量は真の値の近くに分布している。

```R
# 尤度関数(事後分布)
# 正規分布 N(x|a,s)

a0 = 0 # 平均値
s0 = 1 # 標準偏差
N = 500

my_post_norm = function(a,s){
  tmp = sapply(1:length(xn),function(i){
    dnorm(xn[i],mean=a,sd=s,log=T)
  })
  lp = rowSums(tmp)
  exp(lp-max(lp))
}

grid_a = seq(from=-1,to=1,by=0.02)
grid_s = seq(from=0.1,to=2,by=0.02)

set.seed(3)
par(mfrow=c(1,1))
par(mar=c(4,4,2,2))
for(iter in 1:50){
  xn = rnorm(N,mean=a0,sd=s0)
  am = mean(xn)
  sm = sqrt(sum((xn-am)^2)/N)
  z = outer(grid_a,grid_s,my_post_norm)
  msg = sprintf("N(x|a,s), N=%d, iter=%d",N,iter)
  contour(grid_a,grid_s,z,xlab="mean(a)",ylab="sd(s)",main=msg,drawlabels=F)
  points(x=a0,y=s0,type="p",pch=19)
  Sys.sleep(0.1)
}
```

![aic1.png](https://qiita-image-store.s3.amazonaws.com/0/192430/4842fc7d-04c1-c62f-ccb2-4fa11c758e0e.png)

#混合正規分布の尤度関数

- 最尤推定量の変動が激しい。

```R
# 尤度関数(事後分布)
# 混合正規分布 (1-a)*N(x|0,1)+a*N(x|b,1)

a0 = 0.5 # 混合比率
b0 = 0.3 # 他方の平均値
N = 500

my_post_mix = function(a,b){
  tmp = sapply(1:length(xn),function(i){
    log((1-a)*dnorm(xn[i],mean=0,sd=1) +
      a*dnorm(xn[i],mean=b,sd=1))
  })
  lp = rowSums(tmp)
  exp(lp-max(lp))
}

grid_a = seq(from=0,to=1,by=0.02)
grid_b = seq(from=0,to=1,by=0.02)

set.seed(3)
par(mfrow=c(1,1))
par(mar=c(4,4,2,2))
for(iter in 1:50){
  mn = ifelse(runif(N)>a0,0,b0)
  xn = rnorm(N,mean=mn,sd=1)
  z = outer(grid_a,grid_b,my_post_mix)
  msg = sprintf("(1-a)*N(x|0,1)+a*N(x|b,1), N=%d, iter=%d",N,iter)
  contour(grid_a,grid_b,z,xlab="ratio(a)",ylab="mean(b)",main=msg,drawlabels=F)
  points(x=a0,y=b0,type="p",pch=19)
  Sys.sleep(0.1)
}
```

![aic2.png](https://qiita-image-store.s3.amazonaws.com/0/192430/13fd15c6-1b7e-a954-8a0a-c1d82ed384ac.png)

#回帰分析における汎化損失と学習損失の比較

- 予測モデルを$p(x,y|a)=N(y|ax,s_0)U(x)$とする。パラメータ数$d=1$
- $U(x)$は区間[0,1]の一様分布。$s_0$はハイパーパラメータ。$s_0=\sigma$
- 真のモデルを$q(x,y)=p(x,y|a_0)=N(y|a_0x,s_0)U(x)$とする。

```math
\begin{align}
G(a)&=-\iint q(x,y)\log p(x,y|a)dxdy\quad 汎化損失\\
T(a)&=-\frac{1}{N}\sum_{i=1}^N \log p(x_i,y_i|a)\quad 学習損失\\
G(a)&=-\iint N(y|a_0x,s_0)U(x)\log N(y|ax,s_0)U(x)dxdy\\
&=-\iint N(y|a_0x,s_0)\log N(y|ax,s_0)dxdy\\
&=\iint N(y|a_0x,s_0)\left[\frac{1}{2}\log (2\pi s_0^2)+\frac{(y-ax)^2}{2s_0^2}\right]dxdy\\
&=\int\left[\frac{1}{2}\log (2\pi s_0^2)+\frac{(a_0-a)^2x^2}{2s_0^2}+\frac{1}{2}\right]dx\\
&=\frac{1}{2}\log (2\pi s_0^2)+\frac{1}{2}+\frac{(a-a_0)^2}{6s_0^2}\\
T(a)&=\frac{1}{N}\sum_{i=1}^N\left[\frac{1}{2}\log(2\pi s_0^2)+\frac{(y_i-ax_i)^2}{2s_0^2}\right]\\
&=\frac{1}{2}\log(2\pi s_0^2)+\frac{1}{2s_0^2N}\sum_{i=1}^N(y_i-ax_i)^2\\
\hat{a}&=\frac{\sum_{i=1}^N y_ix_i}{\sum_{i=1}^N x_i^2}\quad 最尤推定量\\
E[G(\hat{a})-T(\hat{a})]&=\frac{d}{N}
\end{align}
```

- 最尤推定量の分散

```math
\begin{align}
p(x,y|a)&=N(y|ax,s_0)U(x)\\
\log p(x,y|a)&=\log U(x)-\frac{1}{2}\log(2\pi s_0^2)-\frac{1}{2s_0^2}(y-ax)^2\\
I_{aa}&=-\int\left[
\frac{\partial^2}{\partial a^2}\log p(x,y|a)
\right]p(x,y|a)dxdy\\
&=\int\frac{x^2}{s_0^2}N(y|ax,s_0)U(x)dxdy\\
&=\int_0^1\frac{x^2}{s_0^2}dx=\frac{1}{3s_0^2}\\
\sigma_a^2&=(nI_{aa})^{-1}=\frac{3s_0^2}{n}
\end{align}
```

```R
# 汎化損失と学習損失の差異 N(y|ax,s0)U(x)

# General Loss
my_general_loss = function(a){
  0.5*log(2*pi*s0^2)+0.5+(a-a0)^2/6/s0^2
}

# Training Loss
my_training_loss = function(a){
  sapply(a,function(a1){
    0.5*log(2*pi*s0^2)+mean((yn-a1*xn)^2)/2/s0^2
  })
}

N = 20
a0 = 1
s0 = 0.2
d = 1

set.seed(3)
par(mfrow=c(2,2))
par(mar=c(4,4,2,2))
M = 120
amv = numeric()
amf = numeric()
for(iter in 1:M){
  xn = runif(N)
  yn = a0*xn + rnorm(N,mean=0,sd=s0)

  ymin = min(yn,-0.5)
  ymax = max(yn,1.5)
  msg = sprintf("N(y|ax,s0)U(x), N=%d, iter=%d",N,iter)
  plot(x=xn,y=yn,type="p",xlim=c(0,1),ylim=c(ymin,ymax),main=msg,xlab="X",ylab="Y")
  plot(my_general_loss,0.5,1.5,col="red",ylim=c(-0.7,1),xlab="a",ylab="Loss",
    main="General(red),Training(blue)")
  plot(my_training_loss,0.5,1.5,col="blue",add=T)
  am = sum(yn*xn)/sum(xn^2)
  tm = my_training_loss(am)
  points(x=am,y=tm,type="p",col="blue",pch=1)
  gm = my_general_loss(am)
  points(x=am,y=gm,type="p",col="red",pch=1)

  amv[iter] = am
  amf[iter] = gm -tm
  hist(amv,breaks=seq(from=0.7,to=1.3,by=0.025),ylim=c(0,30),xlab="",main="MLE of a")
  plot(amf,type="l",xlim=c(1,M),xlab="",ylab="General-Training")
  abline(h=1/N,col="red")
  dev.flush()
  Sys.sleep(0.1)
}
```

![aic3.png](https://qiita-image-store.s3.amazonaws.com/0/192430/3ffc2bdd-dc13-1e71-c70c-6fa4c815d286.png)

```R
> c(d/N,mean(amf))
[1] 0.05000000 0.03618879
> c(a0,mean(amv))
[1] 1.000000 1.001853
> c(sqrt(3*s0^2/N),sd(amv))
[1] 0.07745967 0.07874688
> par(mfrow=c(1,1))
> par(mar=c(4,4,2,2))
> qqnorm(amv)
> qqline(amv)
```

![aic7.png](https://qiita-image-store.s3.amazonaws.com/0/192430/802fc948-f656-c09f-4cca-1c3def61675d.png)