R
Stan
ベイズ統計
汎化誤差
ハミルトニアンモンテカルロ法

ベイズ統計「RとStanでベイズ推定の汎化誤差を計算し、理論値と比較する」

はじめに

[追記]はじめて実対数閾値という言葉を聞く人でも分かるような説明(研究分野の紹介(PDF))も作成しました。

今回は、こちらの記事の問題で求めた汎化誤差の理論値にたいして、
RとStanを用いて同じ問題の汎化誤差の実験推定値を計算し、理論値との比較実験を行います。

こちらがRのコードStanのコードです。

予測分布の作り方や汎化誤差については、この本と、その解説をご覧ください。

以下のように問題を設定しました。

問題

データを $ x,y,z\in\mathbb{R} $、パラメータを $ a,b\in\mathbb{R} $とする。

学習モデルを
$$ p(z|x,y,a,b)=\frac{1}{\sqrt{2\pi}}exp(-\frac{1}{2}(z-a^2x-aby)^2) $$
事前分布を
$$ \varphi(a,b)=\begin{cases}
1/4 & |a|,|b|\le1 \\
0 &その他
\end{cases} $$
真の分布を
$$ N(z)N(x)N(y) $$
(ただし、$ N(x) $は標準正規分布)とする。

サンプルがそれぞれn=100,200,400,600,800,1000個与えられたとき、
RとStanでシミュレーションを行い、汎化誤差の実験推定値を計算し、
理論値

$$ \frac{1}{2n}-\frac{1}{nlogn}+o_p(\frac{1}{nlogn}) $$

と比較せよ。
ただし、nがそれぞれの場合について、違うデータセットを適切な数だけ用意し、汎化誤差の実験推定値もその数だけ計算し、
それらをさらにデータセットの出方によって平均したものを最終的な実験推定値とする。
(汎化誤差の実験値は、データセットの出方によってばらつきがあるためである。)

作成したプログラムの説明

まずStanコードの説明をする。

データブロックとパラメータブロックは次のようになる。

BayesEx1
data{
  int<lower=1> N;
  real x[N];
  real y[N];
  real z[N];
}

parameters{
  real a;
  real b;
}

モデルブロックは以下のようになる。

BayesEx1
model{
  a ~ uniform(-1.0, 1.0);
  b ~ uniform(-1.0, 1.0);
  for (n in 1:N) {
    z[n] ~ normal(a*a*x[n]+a*b*y[n], 1.0);
  }
}

モデルブロックのa,bの事前分布は、問題設定より独立であることが分かる。
そのためこのように分けて書くことができる。


次にRコードの説明をしていく。

まずStanを動かして、NUTSハミルトニアンモンテカルロ法によって、事後分布のサンプルを発生させる。

run-BayesEx1
library(rstan)

M = 100 #平均するデータセットの数
N = 100 #サンプルの数
warmup = 200 #HMCのwarmup数
iter = 4000 #HMCのサンプリング数
GN = numeric(M)

for (h in 1:M) {
  #HMCの実行
  x = c(rnorm(N, 0,1))
  y = c(rnorm(N, 0,1))

  z = c(rnorm(N, 0,1))

  data = list(n=N, x=x, y=y, z=z)

  fit = stan(file='ByesEx1.stan', data=data, chains=1, iter=iter, warmup=warmup)

とりあえず、データセットの出方による平均は100個による平均とした。

次にfitから、実際に出たサンプルを取り出す。

run-BayesEx1
  #HMCでサンプリングした事後分布のサンプルを取り出す
  post_param = extract(fit, permuted = FALSE)
  post_a_samp = post_param[1:(iter-warmup), 1, 1]
  post_b_samp = post_param[1:(iter-warmup), 1, 2]

extract()を使うと実際に出たサンプルを取り出せる。

次に確率モデルを事後分布で平均する。

run-BayesEx1
  #予測分布を作る
  p_x_w = numeric(iter-warmup)
  pred_p = function (zp, xp, yp) {
    p_x_w = dnorm(zp, mean =post_a_samp*post_a_samp*xp+post_a_samp*post_b_samp*yp, sd =1)
    return (mean(p_x_w))
  }

予測分布は確率モデルを事後分布で平均したものである。すなはち
事後分布を $ p(w|X^n) $、確率モデルを $ p(x|w) $ としたとき
$$ p^*(x)=\int p(x|w)p(w|X^n)dw $$
が予測分布となるが、この計算は当然解析的にはできないことが多く、サンプルを用いてモンテカルロ積分することになる。(これがベイズ統計の主な特徴の一つである)
それが上のコードである。

予測分布ができたら、次は汎化誤差を作る。

run-BayesEx1
#汎化誤差を作る
  IntN = 5000 #モンテカルロ積分に使うサンプリング数
  G = numeric(IntN)
  Intx = c(rnorm(IntN, 0,1))
  Inty = c(rnorm(IntN, 0,1))
  Intz = c(rnorm(IntN, 0,1))

  for (j in 1:IntN) {
    log_pred_p = log(pred_p(Intz[j], Intx[j], inty[j]))
    log_q = log(dnorm(Intz[j], mean =0, sd =1))
    log = log_q-log_pred_p

    G[j] = log
  }

  GN[h] = mean(G) #汎化誤差

}

汎化誤差を作る際、サンプルが出ている分布で平均する必要があるが、これもモンテカルロ積分で行う。

最後にこれらをデータセットの出方の分だけ繰り返したものを平均する。

run-BayesEx1
EGN11 = mean(GN)

これが汎化誤差の実験推定値となる。

実行結果

実行結果を示す。

Rplot-run-ByesEx1.png

赤:データセットの出方で平均した汎化誤差
黒:データセットの出方で平均する前の100個の汎化誤差
青:理論値(こちらの記事の最後のもの、2つあるうち上側は$ o_p $項を含めたもの、下側は捨てたもの)
である。

数字としては
n=100
理論値:0.002828528(0.005)
実験値:0.002684248

n=200
理論値:0.001556304(0.0025)
実験値:0.001255377

n=400
理論値:0.0008327397(0.00125)
実験値:0.0007161163

n=600
理論値:0.0005727917(0.0008333333)
実験値:0.0005511696

n=800
理論値:0.0004380033(0.000625)
実験値:0.0004297247

n=1000
理論値:0.0003552352(0.0005)
実験値:0.0003289795
(()内は$ o_p $項を含めたもの)

となった。

考察

だいたい理論値と実験値が一致しているように見えるので、ひとまずうまくいっているようである。

しかし、実は何度もこのコードを実行すると、最終的な汎化誤差の実験値(データセットの出方で平均したもの)にばらつきがあることが確かめられる。
そのため、よりよい結果にするためには、データセットの出方による平均をさらに増やす必要があるようである。

また、気づきにくいことだが、
データセットの出方で平均する前の汎化誤差(コード中では変数GN)の中身を見てみると、一部符号がマイナスになっている。
汎化誤差は非負であるはずなのにこうなってしまうのは、平均操作を解析的にではなく、モンテカルロ積分したことによるばらつきであると思われる。

モンテカルロ積分に使うサンプル数を増減させると、それに応じて汎化誤差がマイナスになる割合も変化し、
また、データセットの出方による平均は理論値と重なることから、モンテカルロ積分が原因であると判断したが、気持ち悪いところである。

最後に

今回の例については、RとStanを使ったシミュレーションは、理論値に非常に近い精度を出せることが分かりました。
ベイズ統計の強みの一つに、確率モデル、事前分布、真のモデルによらず、一般的な条件のもとで、汎化誤差の漸近挙動が分かるという点があります。
そのため、さらに多くの例を使って実験してみると面白いかなと思います。

質問などありましたら、気軽にお尋ねください。以上です。