2014-03-05
ぱらぱらめくる『Gaussian estimation: Sequence and wavelet models』

ぱらぱらめくるシリーズ分布推定 Wavelet変換ノンパラメトリック Gaussian sequence model James Stein
こちらで一度ぱらぱらめくったのだけれど、消化不良だったので(それと、色々な意味で必要なので)、もう一度、ぱらぱらめくってみる
PDFはこちら
書きかけ
# わかりたい私のためのGaussian Sequence Model

[ぱらぱらめくるシリーズ][分布推定][Wavelet変換][ノンパラメトリック][Gaussian sequence model][James Stein]ぱらぱらめくる『Gaussian estimation: Sequence and wavelet models』

http://d.hatena.ne.jp/ryamada22/20131127
で一度ぱらぱらめくったのだけれど、消化不良だったので(それと、色々な意味で必要なので)、もう一度、ぱらぱらめくってみる

PDFは
http://statweb.stanford.edu/~imj/GE12-27-11.pdf


基本となる対象はGaussian Sequence Model

n=1個以上の対象があって、その値を$\theta_i;i=1,2,...$とする

これを観察するときに、正規乱数誤差が入るが、この正規乱数は平均0ですべてのiについて等分散であるとする

その観測値を
$$
\begin{equation}
y_i = \theta_i + \epsilon z_i\\
z \sim N(0,1)
\end{equation}
$$
と書ける

これを素直にRでシミュレーションすれば
```{r}
n <- 100
theta <- runif(n)*1 # これはなんでもよい
epsilon <- 2
y <- theta + epsilon * rnorm(n)
```

今、このn個の値の組をn次元変数とすれば、各成分の平均が$\theta_i$で与えられた、全軸に関して等分散なn次元正規乱数ともみなせる。そのつもりでシミュレーションすれば、以下のように、n次元ベクトルを1個発生させることと同じ
```{r}
library(MASS)
# n個の値の観測に入りこむ乱数が相互に独立であることを、mvnorm()の分散共分散行列を対角行列にすることが示している
y <- mvrnorm(1,theta,diag(rep(epsilon,n)))
```

観察データを使って推定しよう

一番単純なのは、観測値そのものを推定値とする方法。最尤推定である。$\hat{\theta_i} = y_i$
この場合のMean Squared Error(MSE)は
$$
\begin{equation}
\sum_{i=1}^n (y_i-\theta_i)^2
\end{equation}
$$

```{r}
sum((y-theta)^2)
```

この値の期待値は
$$
\begin{equation}
n\epsilon^2
\end{equation}
$$
になる
```{r}
n.iter <- 10000
mses <- rep(0,n.iter)
for(i in 1:n.iter){
 tmp.y <- theta + epsilon * rnorm(n)
 tmp.theta <- tmp.y # 推定値はそのまま
 mses[i] <- sum((tmp.y-theta)^2)
}
mean(mses)/n
epsilon^2
```


これとは別の２つの推定法がある

一つ目は少し小さくする方法。
Shrinkage法。

```{r}
s <- runif(1)
n.iter <- 10000
mses <- rep(0,n.iter)
for(i in 1:n.iter){
 tmp.y <- theta + epsilon * rnorm(n)
 tmp.theta <- tmp.y*s # 推定値はsで縮める
 mses[i] <- sum((tmp.theta-theta)^2)

}
mean(mses)/n
epsilon^2
```

観測値より少し小さくしただけなのに、MSEが小さくなりました。

この方法では、$\theta$の値が0近辺のときとそうでないときとで、MSEが小さくできるかどうかが変わってくるので、その様子を確認してみる。

```{r}
mean.theta <- seq(from=-10,to=10,length=50)
means <- rep(0,length(mean.theta))
for(j in 1:length(mean.theta)){
  mses <- rep(0,n.iter)
  for(i in 1:n.iter){
   tmp.y <- theta + mean.theta[j] + epsilon * rnorm(n)
   tmp.theta <- tmp.y*s # 推定値はsで縮める
   mses[i] <- sum((tmp.theta-(theta+mean.theta[j]))^2)
  }
  means[j] <- mean(mses)  
}
plot(mean.theta,means)
abline(h=epsilon^2*n,col=2)
```


このように$\theta$の値の絶対値が小さい範囲では、Shrinkage法の方がMSEが小さくなるが、絶対値が大きくなると返ってMSEは大きくなる様子が見て取れる。

では、絶対値に応じて縮め方を調節するとどうなるだろうか。

以下のように、n個の観測値をn次元ベクトルとみなし、その長さの二乗と、観測誤差の分散とを用いた係数でShrinkageの係数を決めると、n次元ベクトルの長さが0に近い当たりでは、$\epsilon^2 *2$くらいにまでMSEが小さくなり、n次元ベクトルの長さが長くなっても、観測値そのものを推定値としたときのMSEより大きくなることはない、という挙動をすることがわかる。
その係数を式で表すと

$$
\begin{equation}
\frac{(n-2)\epsilon^2}{|y|^2}
\end{equation}
$$
となる。

```{r}
mean.theta <- seq(from=-10,to=10,length=50)
means <- rep(0,length(mean.theta))
for(j in 1:length(mean.theta)){
  mses <- rep(0,n.iter)
  for(i in 1:n.iter){
   tmp.y <- theta + mean.theta[j] + epsilon * rnorm(n)
   s2 <- 1-(n-2)*epsilon^2/sum(tmp.y^2)
   tmp.theta <- tmp.y*s2 # 縮め方を値で調節する
   mses[i] <- sum((tmp.theta-(theta+mean.theta[j]))^2)
  }
  means[j] <- mean(mses)  
}
plot(mean.theta,means,ylim = c(0,epsilon^2*n*1.3))
abline(h=epsilon^2*n,col=2)
abline(h=epsilon^2*2,col=3)
```