Stein's unbiased risk estimator略してSUREという枠組みについてメモしておく。

Steinの補題

標準正規分布に従う確率変数\(Z\sim\mathcal{N}(0,1)\)とする。関数\(f:\mathbb{R}\to\mathbb{R}\)が絶対連続とすると、、

\[ \mathbb{E}[Zf(Z)]=\mathbb{E}[f'(Z)]\]

が成立する。ただし導関数\(f'\)の存在とそれを用いた右辺が発散しないこと、\(\mathbb{E}[f'(Z)]\lt\infty\)は仮定する。これは標準正規分布の確率密度関数を\(\phi(z)=(1/\sqrt{2\pi})\exp(-z^2/2)\)とおいて部分積分を行うことによって次のようにして簡単に証明できる。

\[ \begin{align}\mathbb{E}[f'(Z)]&=\int_{-\infty}^\infty f'(z)\phi(z)dz\\&=\left[f(z)\phi(z)\right]_{-\infty}^\infty+\int_{-\infty}^\infty zf(z)\phi(z)dz\\&=\int_{-\infty}^\infty zf(z)\phi(z)dz=\mathbb{E}[Zf(Z)]\end{align}\]

途中\(\phi'(z)=-z\phi(z)\)という関係を使った。これは\((1/\sqrt{2\pi})\exp(-z^2/2)\)を直接微分してみれば簡単にわかる。

これを\(X\sim\mathcal{N}(\mu,\sigma^2)\)という標準に限らない正規分布に従う確率変数の場合に拡張するには、代わりに\(\psi(x)=(1/\sqrt{2\pi\sigma^2})\exp(-(x-\mu)^2/2\sigma^2)\)という密度関数を使って同様に計算し、最後に\(\psi'(x)=((x-\mu)/\sigma^2)\psi(x)\)という関係式を用いればよく、

\[ \frac1{\sigma^2}\mathbb{E}[(X-\mu)f(X)]=\mathbb{E}[f'(X)]\]

がすぐにわかる。

次にこの多変量正規分布の場合を考える。\(Z^n=(Z_1,\dots,Z_n)\)として\(Z^n\sim\mathcal{N}(0,I_n)\)であるとする。\(I_n\)は\(n\times n\)の単位行列を示している。共分散行列が単位行列であるということは、\(Z^n\)が単に標準正規分布に従う確率変数を\(n\)個並べたものであることを意味している。（\(n\)乗ではないので注意！）

関数\(f:\mathbb{R}^n\to\mathbb{R}\)を考えて、今度は各変数ごとに同じ仮定が成り立つとする。つまり、\(f(z_1,\dots,z_n)\)という関数において、\(z_i\)以外を固定して\(g(z_i)\equiv f(z_1,\dots,z_n)\)と見なしたときの\(g:\mathbb{R}\to\mathbb{R}\)に対して単変量のときと同じ仮定（絶対連続で\(g'\)が存在、さらに\(\mathbb{E}[g'(Z)]\lt\infty\)）が\(i=1,\dots,n\)について成立するとする。

\[ \nabla f=\left(\frac{\partial f}{\partial z_1},\dots,\frac{\partial f}{\partial z_n}\right)\]

と書くようにすると、

\[ \mathbb{E}[Z^nf(Z^n)]=\mathbb{E}[\nabla f(Z^n)]\]

が成立する。Steinの方法（らしい方法）でこれを証明すると、

\[ \begin{align}\mathbb{E}_{z_i}\left[\frac{\partial f(Z^n)}{\partial z_i}\right]&=\int_{-\infty}^\infty\frac{\partial g(z_i)}{\partial z_i}\phi(z_i)dz_i\\&=\int_0^\infty\frac{\partial g(z_i)}{\partial z_i}\phi(z_i)dz_i+\int_{-\infty}^0\frac{\partial g(z_i)}{\partial z_i}\phi(z_i)dz_i\\&=\int_0^\infty\frac{\partial g(z_i)}{\partial z_i}\left(\int_{z_i}^\infty t\phi(t)dt\right)dz_i-\int_{-\infty}^0\frac{\partial g(z_i)}{\partial z_i}\left(\int_{-\infty}^{z_i} t\phi(t)dt\right)dz_i\\&=\int_0^\infty t\phi(t)\left(\int_0^t \frac{\partial g(z_i)}{\partial z_i}dz_i\right)dt-\int_{-\infty}^0 t\phi(t)\left(\int_t^0 \frac{\partial g(z_i)}{\partial z_i}dz_i\right)dt\\&=\int_0^\infty t\phi(t)\left(g(t)-g(0)\right)dt-\int_{-\infty}^0 t\phi(t)\left(g(0)-g(t)\right)dt\\&=\int_{-\infty}^\infty t\phi(t)g(t)dt=\mathbb{E}_{z_i}[z_if(Z^n)]\end{align}\]

途中フビニの定理で積分順序を入れ替えている以外は普通の変形。これは\(\mathbb{E}_{z_i}\)というように\(z_i\)についての期待値しか取っていないので、その他の変数については依然として確率変数。これらについての期待値を取ると、この関係式のどこにも影響なく\(\mathbb{E}_{z_i}\to\mathbb{E}\)とできる。よって

\[ \mathbb{E}\left[\frac{\partial f(Z^n)}{\partial z_i}\right]=\mathbb{E}[z_if(Z^n)]\quad(i=1,\dots,n)\]

がわかった。これは最初の式のベクトル表記を\(i=1,\dots,n\)ごとに書き下したものになっている。

一変数の際と同様にこれを標準でない場合\(X^n\sim\mathcal{N}(\mu,\sigma^2I_n)\)にも一般化すると以下の式になることがわかる。

\[ \frac1{\sigma^2}\mathbb{E}[(X^n-\mu)f(X^n)]=\mathbb{E}[\nabla f(X^n)]\]

SUREで用いる形式に書き換えるために、\(\mathbb{E}[f(X^n)]=0\)を仮定して、\(\mu=\mathbb{E}[X^n]\)を用いると

\[ \begin{align}\mathbb{E}[(X^n-\mu)f(X^n)]&=\mathbb{E}[\;(X^n-\mathbb{E}[X^n])(f(X^n)-\mathbb{E}[f(X^n)])\;]\\&=\mathrm{Cov}[\; X^n,f(X^n)\;]\end{align}\]

より

\[ \frac1{\sigma^2}\mathrm{Cov}[\;X^n,f(X^n)\;]=\mathbb{E}[\nabla f(X^n)]\]

がわかる。ここでの\(\mathrm{Cov}\)はベクトルごとの意味であり、書き下すと

\[ \frac1{\sigma^2}\mathrm{Cov}[x_i, f(X^n)]=\mathbb{E}\left[\frac{\partial f(X^n)}{\partial x_i}\right]\quad(i=1,\dots,n)\]

の意味であることに注意が必要。この補題の結果を用いることによってSURE：Steinの不偏リスク推定というものを考えることが出来る。

SURE: Stein's unbiased risk estimator

\(n\)個の観測点\(y_1,\dots,y_n\)が次に従って生成されている状況を考える。

\[ y_i=\mu_i+\sigma_n z_i\quad(i=1,\dots,n)\]

\(\mu_i\)が未知で\(\sigma_n\)は既知であるとする。ここで\(z_i\sim\mathcal{N}(0,1)\;\)i.i.d. で言い換えると\(z_i,\;i=1,\dots,n\)は独立に標準正規分布に従う。以上は\(y=(y_1,\dots,y_n)\), \(\mu=(\mu_1,\dots,\mu_n)\)を使って単に\(y\sim\mathcal{N}(\mu,\sigma_n^2 I_n)\)とも書ける。

観測\(y\)から未知の\(\mu\)を推定する手続きを関数\(\widehat\mu:\mathbb{R}^n\to\mathbb{R}^n\)を使って表す。\(\widehat\mu(y)\)は\(y\)をもとに\(\mu\)を推定した結果を表す。このとき、以下に示すリスク関数を考える。

\[ R=\mathbb{E}[(\mu-\widehat\mu(y))^2]\]

ベクトル\(y\)の2乗は\(y^2=y^\top y\)によって定義する。すると、これは未知である真の\(\mu\)と観測\(y\)から推定によって得た\(\widehat\mu(y)\)の2乗誤差になる。\(z_1,\dots,z_n\)に関する期待値を取っている（他にランダムな変数は無いため、表記上は省略されている）ため、毎回の推定におけるその量の良さではなく推定手続きそのものの良さを測っているといえる。

今なんらかの方法でリスクを推定したい。単に期待値を外すだけならば\(Z^n\)のindependent copy（乱数を振り直した複製）を\(M\)個用意してその観測を\(y^{(m)},\;m=1,\dots,M\)とし

\[ \widehat{R}=\frac1M\sum_{m=1}^M(\mu-\widehat\mu(y^{(m)}))^2\]

などとすればよいと思うかもしれないが、\(\mu\)は未知なのでそもそもindependent copyを作れないし\(\mu\)を含むこの量の計算も出来ない（モチベーションとしてはありえる）。そこで\(R\)を次のように変形する。以下\(\widehat\mu(y)=\widehat\mu\)として\(y\)の依存性は省略する。

\[ \begin{align}\mathbb{E}[(\mu-\widehat\mu)^2]&=\mathbb{E}[(\mu-y+y-\widehat\mu)^2]\\&=\mathbb{E}[(\mu-y)^2+2(\mu-y)^\top(y-\widehat\mu)+(y-\widehat\mu)^2]\\&=\mathbb{E}[(\mu-y)^2]+2\mathbb{E}[(\mu-y)^\top(y-\widehat\mu)]+\mathbb{E}[(y-\widehat\mu)^2]\\&=n\sigma_n^2+\mathbb{E}[(y-\widehat\mu)^2]+2\mathbb{E}[(\mu-y)^\top(y-\mu+\mu-\widehat\mu)]\\&=n\sigma_n^2+\mathbb{E}[(y-\widehat\mu)^2]-2\mathbb{E}[(\mu-y)^2]+2\mathbb{E}[(y-\mu)^\top(\widehat\mu-\mu)]\\&=-n\sigma_n^2+\mathbb{E}[(y-\widehat\mu)^2]+2\sum_{i=1}^n\mathrm{Cov}[y_i,\widehat\mu_i]\end{align}\]

と展開できる。ただし\(\mathbb{E}[(\mu-y)^2]=\mathbb{E}[(\sigma_nz)^2]=n\sigma_n^2\)や、\(\mathbb{E}[y]=\mu\)が仮定から従うことを使っている。さらに、ここではわかりやすさのため\(\mathbb{E}[\widehat\mu(y)]=\mu\)（推定量の不偏性）も仮定している。

第3項は一見して計算可能だがこの計算には\(\mu\)が必要になる。ところがこの部分はSteinの補題のところで最後に示した結果をそのまま用いると\(\mu\)を知ることなく計算できる。これを踏まえて期待値を外した以下のリスク推定量

\[ \widehat R=-n\sigma_n^2+(y-\widehat\mu)^2+2\sigma_n^2\sum_{i=1}^n\frac{\partial\widehat\mu_i(y)}{\partial y_i}\]

をStein's unbiased risk estimatorという。略してSUREといわれestimatorではなくestimationの意味でも使われる（推定量と推定の違い）。実用的には\(n\)の大きさが変化した時にリスクの大きさが変わってしまうのは都合が悪いため、最初から\(R=\mathbb{E}[(\mu-\widehat\mu)^2]/n\)を考えていたとすると、SUREは

\[ \widehat{R}=-\sigma_n^2+\frac1n\sum_{i=1}^n(y_i-\widehat\mu_i)^2+\frac{2\sigma_n^2}{n}\sum_{i=1}^n\frac{\partial\widehat\mu_i(y)}{\partial y_i}\]

という形になり、既知の量だけから計算することができる。

第2項は推定誤差になっており、文脈によっては訓練誤差ともいわれる。極端な話\(\widehat\mu(y)=y\)というルールを用いればこの部分はゼロに出来るが、その場合は第3項の偏微分の箇所は全て1になり、結果\(\widehat{R}=\sigma_n^2\)という結果になる。これに関してはSteinのパラドックスという面白い話が知られている。

第2項と第3項のバランスによってもっと小さい\(\widehat{R}\)が実現出来たりできなかったりするわけである。そういう訳でSUREはモデル選択に用いることが出来て、この場合後半の第3項はモデルの自由度を定量化しているとみるのが自然である。実際、

\[ \mathrm{df}(\widehat\mu)=\frac1{\sigma_n^2}\sum_{i=1}^n\mathrm{Cov}[\;y_i,\widehat\mu_i(y)\;]\approx\sum_{i=1}^n\left[\frac{\partial \widehat\mu_i(y)}{\partial y_i}\right]\]

という量でもって"Degrees of freedom"、モデルの自由度\(\mathrm{df}(\widehat\mu)\)を定義することがありSUREではこれを最右辺の量でもって推定する。

以下違うページに続く。

参考

Ryan Tibshirani (with Larry Wasserman), Stein's Unbiased Risk Estimate , Statistical Machine Learning, Spring 2015, (http://www.stat.cmu.edu/~larry/=sml/stein.pdf)