スタイン推定量を確かめてみる - 統計コンサルの議事メモ

いきなりですが、最近購入したベイズモデリングの世界を読んでいて非常に面白い話題を発見しました。

作者: 伊庭幸人
出版社/メーカー: 岩波書店
発売日: 2018/01/18
メディア: 単行本（ソフトカバー）
この商品を含むブログ (4件) を見る

P118から始まるスタイン推定量についてなのですが、その直感的な意味について本文から引用すると：

それぞれの $y_i$ を、全部の{ $y_i$ }の平均値 $\frac{1}{n}\Sigma_{i=1}^{n}y_i$ の方向に $a$ だけ引っ張ってやる

ことで、 $y_i$ そのものを推定量とするよりも、パラメータの真値 $\hat{\theta_{i}}$ との二乗誤差の期待値が小さくなるそうです。これはいわゆる縮小推定量の一種で、引っ張りの程度 $a$ をデータから適応的に求めるのがこの推定量のキモなのだそうです。

スタイン推定量 $\hat{\theta_{i}^{S}}$ の定義自体はそれほど難しいものでもなく、以下のような式によって表されます：

$\displaystyle \hat{\theta_{i}^{S}} = (1-a)y_{i} + \frac{n}{a}\Sigma_{i=1}^{n}y_i$

ここで、

$\displaystyle a = \frac{\sigma^{2}}{s^{2}}, \ s^{2} = \frac{1}{n-3}\Sigma_{i=1}^{n}(y_i - \frac{1}{n}\Sigma_{j=1}^{n}y_i)^{2}$

です。この記事では、スタイン推定量が本当に $y_i$ よりも二乗誤差が小さくなるのか、以下のように実験してみます。

実験内容

内容としては非常に簡単なもので、平均と分散が共に既知である正規分布からいくつかのサンプルを生成します。そのサンプルをそのまま用いた場合と、スタイン推定量を用いた場合とで、真のパラメータとの二乗誤差の期待値：

$\mathbb{E}[\Sigma_{i=1}^{n}(\theta_{i}^{*}({y_i}) - \theta_i)^{2}$ ]

にどれほど差が生じるかを確認しています。

パラメータ設定

実験では、一回の試行につき平均が{1, 2, ..., n}、分散が1である正規分布に従うサンプルをn個生成することとしました。このn個のデータから各サンプルのスタイン推定量を求め、真値である{1, 2, ..., n}との二乗誤差を計算します。同時にデータそのものを用いた場合の二乗誤差も計算しておきます。上記の実験を1000回繰り返し、最終的にどちらが小さい値となるかを比較しました。

n   <- 10
mu  <- 1:n
sig <- 1

K   <- 1000
res <- matrix(0, K, 2)

実験開始

上記のパラメータに基づき、以下のように計算を行います。

set.seed(123)
for (i in 1:K) {
   ## サンプルを発生させる
   tmp <- rnorm(n, mu, sig)
   
   ## スタイン推定値を求める
   s2  <- (1/(n-3)) * sum((tmp - mean(tmp))^2)      
   a   <- sig / s2
   s_i <- (1-a)*tmp + (a/n)*sum(tmp)

   ## 結果を保存する
   res[i, 1] <- sum((s_i - mu)^2)/n # スタイン推定値を用いた二乗誤差の期待値
   res[i, 2] <- sum((tmp - mu)^2)/n # データそのものを用いた二条誤差の期待値
}

結果

それでは結果です。

> mean(res[, 1] / res[, 2])
[1] 0.9618857

> sum(res[, 1] > res[, 2])
[1] 310

スタイン推定量を用いた場合、データそのものの場合と比較して二乗誤差が96%と小さくなっていることが確認できます。

終わりに

今回の実験は以上で終了です。ではこれの何が面白いのかと言うと、スタイン推定量と階層ベイズモデルとの関連性についてなんですね。現在、別に書こうとしている記事で階層ベイズモデルに触れているのですが、その説明として

パラメータを生成する分布を仮定することで、パラメータに緩やかな縛りをかける

といった言い回しを考えています。ここでふと「階層ベイズでは（緩やかな縛りをかけることで）推定値にバイアスを与えているわけだが、それにも関わらずモデルとして良くなるというのは、一体どういう意味なのだろう？」という疑問が浮かんできました。バイアスを与えているわけなのだから、誤差が大きくなってモデルとしては劣化しそうな気がしませんか？

その時にこのスタイン推定量のことを思い出しました。つまり、階層ベイズではパラメータに対する分布を仮定することで推定値にバイアスを与えているが、スタイン推定量の性質を以って全体的には誤差を小さくすることができているということではないかと¹。

書籍にはスタイン推定量が誤差を減少させる証明も付いていますが、まだ読み込めていませんので上記の理解が間違っているかもしれません。しかしこの推定量自体も非常に面白く、紹介したいと思ったので実験してみた次第です。