単純モデルについてWAIC・WBICを算出

#概要

- 単純モデル$(Y\sim N(a,1/s))$について、WAIC、AIC、BIC、Gn、WBIC、Fを算出。
- Gibbs samplingを使用。
- Gn(汎化損失)については数値積分を使用。
- 「渡辺澄夫 [広く使える情報量規準(WAIC)](http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/waic2011.html)」
- 「渡辺澄夫 [広く使えるベイズ情報量規準(WBIC)](http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/wbic2012.html)」

#WAICの計算式

```math
\begin{align}
p(y|a,s)&=\sqrt{\frac{s}{2\pi}}\exp\left(-\frac{s}{2}(y-a)^2\right)\quad モデル\\
\phi(a,s|h,r)&\propto s^r\exp\left(-\frac{hs}{2}(a^2+1)\right)\quad 事前分布(r=0,h=0.001)\\
p(a,s|D)&\propto \phi(a,s|h,r)\prod p(Y_i|a,s)\quad 事後分布\\
&\propto s^{r+\frac{n}{2}}\exp \left(-\frac{sA}{2}(a-B)^2-\frac{C}{2}s\right)\\
&A=n+h,B=\sum Y_i/A,C=Y^tY+h-AB^2\\
Z&=\iint s^{r+\frac{n}{2}}\exp \left(-\frac{sA}{2}(a-B)^2-\frac{C}{2}s\right)da ds\\
&=\sqrt{\frac{2\pi}{A}}\frac{\Gamma(\alpha)}{\beta^\alpha},\alpha=r+\frac{n+1}{2},\beta=\frac{C}{2}\\
p(y)&=E_w[p(y|a,s)]\quad 予測分布\\
q(y)&=p(y|a_0,s_0),a_0=1,s_0=10\quad 真のモデル\\
S&=-E_Y[\log q(y)]=\frac{1}{2}+\frac{1}{2}\log \frac{2\pi}{s_0}\quad 真の分布のエントロピー\\
G_n&=-E_Y[\log p(y)]\quad 汎化損失\\
T_n&=-\frac{1}{n}\sum\log p(Y_i)\quad 学習損失\\
V_n&=\sum\left(E_w[(\log p(Y_i|a,s))^2]-E_w[(\log p(Y_i|a,s))]^2\right)\quad 汎関数分散\\
\mathrm{WAIC}&=T_n+V_n/n
\end{align}
```

#WBICの計算式

```math
\begin{align}
\mathrm{WBIC}&=-\iint\sum\log p(Y_i|a,s)p_\beta(a,s|D)dads\quad(\beta=1/\log n)\\
p_\beta(a,s|D)&\propto \phi(a,s|h,r)\prod p(Y_i|a,s)^\beta\\
&\propto s^{r+\beta\frac{n}{2}}\exp \left(-\frac{sA}{2}(a-B)^2-\frac{C}{2}s\right)\\
&A=\beta n+h,B=\beta\sum Y_i/A,C=\beta Y^tY+h-AB^2\\
F&=-\log\iint\phi(a,s|h,r)\prod p(Y_i|a,s)dads\quad ベイズ自由エネルギー
\end{align}
```

#Rスクリプト

```
# lz = log int s^r exp(-s/2*(Aa^2+C))
my_get_lz <- function(r,A,C){
    log(2*pi/A)/2+lgamma(r+1/2)-log(C/2)*(r+1/2)
}
my_get_F <- function(r,h,y){
    n = length(y)
    a = n+h
    b = sum(y)/a
    c = sum(y*y)+h-a*b*b
    F = my_get_lz(r,h,h) + log(2*pi)*n/2 - my_get_lz(r+n/2,a,c)
}
# Gibbs sampling
my_gibbs_sam <- function(m,r,h,y,beta){
    n = length(y)
    a = beta*n+h
    b = beta*sum(y)/a
    c = beta*sum(y*y)+h-a*b*b
    shape = r+beta*n/2+1
    rate = c/2
    ap = b
    sp = shape/rate
    al = 0
    sl = 0
    for(i in 1:m){
        for(j in 1:10){
            ap = rnorm(1,mean=b,sd=1/sqrt(sp*a))
            rate = a/2*(ap-b)*(ap-b)+c/2
            sp = rgamma(1,shape=shape,rate=rate)
        }
        al[i] = ap
        sl[i] = sp
    }
    list(a=al,s=sl)
}

# waic
tm = 100
t0 <- proc.time()
reta = t(sapply(1:tm,function(seed){
    set.seed(seed)
    n = 100
    a0 = 1
    s0 = 10
    sd0 = 1/sqrt(s0)
    Y = rnorm(n,mean=a0,sd=sd0)
    r = 0
    h = 0.001
    # Gibbs sampling
    m = 10000
    sam = my_gibbs_sam(m,r,h,Y,1)
    al = sam$a
    sl = sam$s
    # 汎化損失
    yl = a0-sd0*6
    yu = a0+sd0*6
    hp = (yu-yl)/100
    Gn = 0
    for(j in 0:100){
        y2 = yl+hp*j
        pf = mean(dnorm(y2,mean=al,sd=1/sqrt(sl)))
        Gn = Gn - log(pf) * dnorm(y2,a0,sd0) * hp
    }
    # 真の分布のエントロピー
    S = 1/2+1/2*log(2*pi/s0)
    # 汎関数分散
    Vn = 0
    vl = vector(length=n)
    for(j in 1:n){
        vl[j] = mean(dnorm(Y[j],mean=al,sd=1/sqrt(sl)))
        v = dnorm(Y[j],mean=al,sd=1/sqrt(sl),log=T)
        Vn = Vn+mean(v*v)-(mean(v))^2
    }
    # 学習損失
    Tn = -mean(log(vl))
    waic = Tn+Vn/n
    # 最尤推定(の近似値)
    am = mean(al)
    sm = mean(sl)
    np = 2
    aic = -sum(dnorm(Y,mean=am,sd=1/sqrt(sm),log=T))*2+np*2
    bic = -sum(dnorm(Y,mean=am,sd=1/sqrt(sm),log=T))*2+np*log(n)
    # output
    c(n=n,S=S,Gn=Gn,Tn=Tn,Vn=Vn,waic=waic,aic=aic,bic=bic)
}))
proc.time() - t0

# wbic
tm = 100
t0 <- proc.time()
retb = t(sapply(1:tm,function(seed){
    set.seed(seed)
    n = 100
    a0 = 1
    s0 = 10
    sd0 = 1/sqrt(s0)
    Y = rnorm(n,mean=a0,sd=sd0)
    r = 0
    h = 0.001
    beta = 1/log(n)
    # Gibbs sampling
    m = 10000
    sam = my_gibbs_sam(m,r,h,Y,beta)
    al = sam$a
    sl = sam$s
    wbic = -sum(sapply(1:n,function(j){mean(dnorm(Y[j],mean=al,sd=1/sqrt(sl),log=T))}))
    f = my_get_F(r,h,Y)
    # output
    c(n=n,wbic=wbic,f=f)
}))
proc.time() - t0
```

#計算結果

```
head(cbind(reta,retb))

       n        S        Gn        Tn       Vn      waic      aic      bic   n     wbic        f
[1,] 100 0.267646 0.2860429 0.1553681 2.067547 0.1760436 35.05913 40.26947 100 20.11540 25.44235
[2,] 100 0.267646 0.2889707 0.4131372 1.674157 0.4298788 86.24508 91.45542 100 45.75631 51.28737
[3,] 100 0.267646 0.2925912 0.1090112 1.643863 0.1254498 25.45273 30.66307 100 15.30152 20.59201
[4,] 100 0.267646 0.2806845 0.1739265 1.770756 0.1916341 38.51656 43.72690 100 21.81383 27.18913
[5,] 100 0.267646 0.2707372 0.2070117 1.926296 0.2262747 45.27221 50.48255 100 25.17268 30.60127
[6,] 100 0.267646 0.2690678 0.2964459 1.859406 0.3150399 63.11569 68.32604 100 34.19786 39.60833
```

#計算結果の確認

- waicの平均値と汎化損失(Gn)の平均値はほぼ同じとなっている。
- waic=aic/(2n)、wbic=bic/2の関係がほぼ成り立っている。
- wbicとf(ベイズ自由エネルギー)はほぼ一致するはずだが、fの方が一貫して5.4くらい大きい。
- 原因が分からない。ベイズ自由エネルギーの定義が異なっている?