数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

『わたしのためのJames-Stein Estimator』

James-Stein Estimator R epub knitr
  • こちらでGaussian estimation: Sequence and wavelet modelsのおさらいをしようとしている
  • その手前で「複数(たくさん)の観察がなされたときに、それらを一括してデータマイニングする」一つの方法としてJames-Stein Estimatorが出てくる
  • Gaussian Estimationでは、James-Stein Estimatorを「全レコードについて、あるルールで縮める手法」として取り上げ、それ以外の手法を「部分的に圧縮する」というように考え、両者を同じ枠組みで考えている。全体を圧縮するか、圧縮に濃淡を入れるか、圧縮にゼロ・イチを入れるかなどの違いはあるが、意味合いは同じ、ということである。「意味のあるところだけを取り出す」手法としてスプラインとかフーリエ変換を介した雑音除去、ウェーブレット変換を介した雑音除去なども、ゼロ・イチ圧縮の手法として併せて論じている
  • というわけで、その一番手前のJames-Stein Estimatorについて
  • ウェブ上PDF(こちら)とそれのモノグラフ化したもの(以下のアマゾン書籍)とがあり、その第１章がJames-Stein Estimatorについてなのだが、「理論的」で、「わかりたいだけの私」にはつらい

Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction (Institute of Mathematical Statistics Monographs)

Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction (Institute of Mathematical Statistics Monographs)

  • というわけで、同書の第１章の内容をR的になぞってみたのが以下。
  • epub化したので、こちらに登場予定だけれど、Rmdファイルなので、ここを参考にフリーでhtml化・epub化は可能です

James-Stein Estimator?理論式よりシミュレーションで? 統計遺伝学

James-Stein Estimator?理論式よりシミュレーションで? 統計遺伝学

# 私のためのJames-Stein Estimator

## はじめに

複数の対象に興味があるときに、複数の対象について観察がなされたとき、その観察結果のすべてを使って対象の推定をすると、個々の対象について個別に推定するよりもよい、という話。

ごく限定的な状況からはじめて、順々にこみいった状況に進めていくことにする。

理論的な式変形は他の本にいくらでもあるから、なるべくそれは使わずに、シミュレーションで納得することを基本とする。

参考文献１： http://statweb.stanford.edu/~omkar/329/chapter1.pdf (Chapter1)

参考文献２： http://en.wikipedia.org/wiki/James%E2%80%93Stein_estimator 

## James-Stein Estimatorの心

対象数をnとする。
それらは、何かしら期待値$M$を持つが、実際には、そこから分散$A$でずれているものとする。
そのような値を「真値」と呼び、$\mathbf{\theta}=(\theta_1,...,\theta_n)$という値になっているとする。

ここで、n個を観察すると、観察誤差$\mathbf{z}=(z_1,...,z_n)$が入る。
観察誤差は期待値が$m=0$で分散が$a$であるとする。

$$
\begin{equation}
y_i = \theta_i + z_i
\end{equation}
$$


このようにして観察された値は$\mathbf{y}=(y_1,...,y_n)$であるとする。

この$\mathbf{y}$の値から、次のようにして真値$\mathbf{\theta}$の推定値$\hat{\mathbf{\theta}}$を算出することとする。

このとき、$\mathbf{y} = (y_1,...)$から、$\mathbf{\theta}=(\theta_1,...)$の値を次のように推定することにして、推定値を$\hat{\mathbf{\theta}}=(\hat{\theta_1},...)$と書くことにする。

$$
\begin{equation}
\hat{\theta_i} = k(y_i - M) + M\\
k = \frac{A}{A+a}
\end{equation}
$$

この心は、

>n 個の対象にそれぞれ色々な値が観測されるけれど、ひとまず、すべては根っこのところでは等しく値がMであって、そこから、「真のずれ」と「観測誤差としてのずれ」の２段階のずれがあると考える。今、観察値の根っこの値からのずれのうち、「真のずれ」の割合だけに減らしてやったら、それが「真のずれだけを含む真値」の推定値としてよいだろう

というものである。


本当は観察値$\mathbf{y}$のみがあって、そこから推定したいわけだけれど、今は、$M$も$A$も$a$の値もわかっているものとする。

このような推定をする関数は
```{r}
js.1 <- function(y,M,A,a){
  n <- length(y)
  k <- A/(A+a)
  ret <- k*(y-M) + M
  return(ret)
}
```

と書ける。


## やはり最初は標準的な正規分布から

対象は１つの値を持ち、それは標準正規分布に従うものとする。
標準正規分布は、平均$M=0$,分散$A=1$なので、
$$
\begin{equation}
\theta_i \sim N(M=0,A=1)\\
i = 1,2,...,n
\end{equation}
$$

n個の値は相互に独立と考えている。

観察すると誤差が入る。
その誤差も標準正規分布とする。

誤差を、その平均$m=0$と分散$a=1$として
$$
\begin{equation}
z_i \sim N(m=0,a=1)
\end{equation}
$$
とすれば、観察値は
$$
\begin{equation}
y_i = \theta_i + z_i
\end{equation}
$$



適当にnを決めて、これを何度もやってみて、確かに、推定値が真値のよい推定になっていることを確かめる。

「よい推定」であるかどうかの基準として、
$$
\begin{equation}
R = \sum_{i=1}^n (\hat{\theta_i}-\theta_i)^2
\end{equation}
$$
を用いることにする。

「よい」ことを示すために、比較する対照として、観察値そのものを推定値とした場合($\hat{\theta_i}_{ML}=y_i$)を用いることにする。
この推定値のことは「最尤推定(ML:Maximum Likelihood)値」と呼ぶ。

```{r}
# 繰り返し試行回数
n.iter <- 1000
# 対象数
n <- 50
# 真値の平均値M、誤差の平均値m
M <- m <- 0
# 真値・誤差の分散
A <- a <- 1
# 試行ごとの、「ずれ」を記録する
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A)) # 分散Aなら標準偏差sqrt(A)
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  # 最尤推定値はyそのもの
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate",ylab="ML estimate")
abline(0,1,col=2)
```

ほとんどの常に、$\hat{\theta}$の場合の方が$\hat{\theta}_{ML}$の場合よりも、ずれの大きさRの値が小さくなっていることがわかる(両者の大きさが等しいのが赤の補助線である)。

## 標準正規分布ではない場合

前節の例は真値が標準正規分布であり、誤差も標準正規分布であるという、特殊な場合であった。

まずは、真値を平均0($M=0$)、分散が0ではない($A \ne 1$)に変えて、同様のことをやってみる。

```{r}
n.iter <- 1000
n <- 50
M <- m <- 0
A <- 3 # 真値の分散を変更
a <- 1
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="ML estimate",ylab="JS estimate")
abline(0,1,col=2)
```

同様の傾向がある。

では、今度は、真値の平均を0でなくしてみる。
分散も1ではないものとする。
```{r}
n.iter <- 1000
n <- 50
M <- 10 # 真値の平均を変更
m <- 0
A <- 3 # 真値の分散を変更
a <- 1
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="ML estimate",ylab="JS estimate")
abline(0,1,col=2)
```


ではさらに、誤差の分散も変えてみる。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
m <- 0
A <- 5 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate",ylab="ML estimate")
abline(0,1,col=2)
```

## 正規分布ではない分布ではどうか

真値の分布を正規分布ではなく、一様分布にしてみる。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均
A <- (M*2)^2/12 # 一様分布の分散
m <- 0
a <- 7 # 誤差の分散を変更
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  # 平均Mの一様分布
  theta <- runif(n)*2*M
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate",ylab="ML estimate")
abline(0,1,col=2)
```

真値の分布を正規分布ではなく、指数分布にしてみる。

```{r}
n.iter <- 1000
n <- 50
M <- 2 # 真値の平均
A <- M^2 # 真値の分散
m <- 0
a <- 7 # 誤差の分散を変更
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  # 平均Mの指数分布
  theta <- rexp(n,1/M)
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate",ylab="ML estimate")
abline(0,1,col=2)
```

では、今度は、誤差項を正規分布から変えてみる。
真値は正規分布に戻そう。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
A <- 5 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  # 分散がaの指数分布乱数
  z <- rexp(n,1/sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate",ylab="ML estimate")
abline(0,1,col=2)
```

最後に、真値も正規分布ではなく、誤差項も正規分布ではない場合をやっておく。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
A <- (M*2)^2/12 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  # 平均Mの一様分布
  theta <- runif(n)*2*M
  # 分散がaの指数分布乱数
  z <- rexp(n,1/sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate",ylab="ML estimate")
abline(0,1,col=2)
```

## 観測値だけから推定したい。

上記の例では、確かによりよい推定ができることがわかった。

しかしながら、推定するにあたり、観察値$\mathbf{y}$の他に、真値の期待値$M$と、真値の$M$からの分散$A$と誤差項の分散$a$とを用いていた。

実際には$M,A,a$のいずれも不明であることがほとんどである。

それらを何かしらで代用しなくては、観察値だけが得られる場合には役に立たない。

## 真値の平均

$M,A,a$のうち、$M$は

$$
\begin{equation}
\hat{M} = \frac{\sum_{i=1}^n y_i}{n}
\end{equation}
$$

で代用するのは悪くない。

実際、$M$を与える代わりに、観測値の平均値$\hat{M}$で代用する関数は以下のようになる。

```{r}
js.2 <- function(y,M=mean(y),A,a){
  n <- length(y)
  k <- A/(A+a)
  ret <- k*(y-M) + M
  return(ret)
}
```

いくつか見てきた例のうちの、真値が平均が0でない正規分布に従い、誤差の分散も1ではない場合で見てみよう。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
m <- 0
A <- 5 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
Rs <- Rs2 <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  theta.hat.2 <- js.2(y,A=A,a=a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs2[i] <- sum((theta-theta.hat.2)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs2,pch=20,asp=1,xlab="Known M",ylab="Estimated M")
abline(0,1,col=2)
```

横軸に$M$がわかっている場合、縦軸に$M$を観察値から推定した場合を取っている。
観測値の平均で代用する場合の方が、少し当てはまりが悪くなるが、ほぼ$y=x$の線に乗っており、代用も悪くないことが見て取れる。

## 真値の分散と誤差の分散との和

### 真値も誤差も正規分布の場合

真値が平均$M=0$、分散$A$で、誤差が平均$m=0$、分散$a$の正規分布に従うとき、観測値は、平均0、分散$A+a$の正規分布に従う。

$$
\begin{equation}
y \sim N(0,A+a)
\end{equation}
$$

このような場合には、真値の分散と誤差の分散の和は、観測値の分散から推定できる。

確かめてみよう。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
m <- 0
A <- 5 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
vars <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  vars[i] <- var(y)
}
mean(vars)
A+a
hist(vars)
abline(v=A+a,col=2)
```

$A+a$の値の推定値の分布と実際の推定値(赤い線)とを示した。


ばらつきはあるが、推定値として使えることがわかる。

$A+a$の期待値の推定はできた。

しかしながら、推定するにあたって用いる補正係数は、

$$
\begin{equation}
\hat{\theta_i} = k(y_i - M) + M\\
k = \frac{A}{A+a}
\end{equation}
$$

に示すように、$\frac{A}{A+a}$であって、$A+a$ が分母に来ている。

この値の推定値を得るには、$\frac{A}{A+a}$自体、もしくは$1-\frac{A}{A+a}=\frac{a}{A+a}$が推定できればよく、$A$または$a$の推定値がわかるのであれば$\frac{1}{A+a}$の推定値がわかれば良いことになる。

一般に、ある分布に従う確率変数$X$の期待値と$\hat{X}$と、その逆数$Y=\frac{1}{X}$の期待値$\hat{Y}$とは$\hat{Y} \ne \frac{1}{\hat{X}}$であるから、$A+a$の期待値がわかるだけでは、目的を達することはできない。

念のため、変数の期待値と変数の逆数の期待値との関係について、簡単な例を示しておく。

```{r}
n <- 1000
X <- runif(n)
mean(X)
mean(1/X)
1/mean(X)
```


従って、$A+a$の推定値が得られるだけではだめで$\frac{1}{A+a}$の推定ができなくてはならないことの確認ができた。

次のように考える。

n次元標準正規分布の原点からの距離の二乗の分布を自由度nの$\chi^2$分布と言う。


n個の対象の観測値$\mathbf{y}$を長さnのベクトルと見たとき、これは、分散が$A+a$の正規分布に従う乱数と見ることができる。
従って、そのベクトルの長さの２乗は自由度nの$\chi$分布を$A+a$倍分だけ$\chi^2$値についてスケール変換した分布に従う。


念のため、確かめておく。

```{r}
n.iter <- 1000
n <- 10
S <- 5
X <- matrix(rnorm(n.iter*n,0,sqrt(S)),ncol=n)
L <- apply(X^2,1,sum)

Y <- rchisq(n.iter,n)
plot(sort(L),sort(Y)*S,xlab="Squared Length of Vector",ylab="Chi-squared dist with df=n")
abline(0,1,col=2)
```

従って、自由度nの$\chi^2$分布に従う確率変数の逆数の期待値がわかればよい。

逆$\chi^2$分布と呼ばれる分布がそれで、その期待値は$\frac{1}{n-2}$となることが知られている。

従って、$\frac{1}{\frac{|\mathbf{y}|^2}{A+a}}=\frac{A+a}{|\mathbf{y}|^2}$は自由度nの逆$\chi^2$分布に従うことから、この期待値は$\frac{1}{n-2}$に従う。

$$
\begin{equation}
E(\frac{A+a}{|\mathbf{y}|^2}) = \frac{1}{n-2}
\end{equation}
$$

これより
$$
\begin{equation}
\frac{1}{A+a} = E(\frac{n-2}{|\mathbf{y}|^2})
\end{equation}
$$

確かめてみよう。
```{r}
n.iter <- 1000
n <- 10 # 要素数
A <- 5
a <- 2
Y <- matrix(rnorm(n.iter*n,0,sqrt(A+a)),ncol=n)
L <- apply(Y^2,1,sum)　# 長さnのベクトルの長さの二乗
Z <- (n-2)/L # この標本の平均が1/(A+a)になっていてほしい
hist(Z)
abline(v=1/(A+a),col=2)
mean(Z)
1/(A+a)
```

試行を繰り返すとばらつくが、真値(赤線)周辺に分布している様子がわかり、実際、全試行での平均値は真値に近い。

## 真値の平均の推定が及ぼす影響

ここで行った$\frac{1}{A+a}$の推定は$\mathbf{y} \sim N(0,A+a)$と、0を中心とした正規分布を仮定していた。
0を中心とした正規分布の二乗の値の分布である$\chi^2$分布が使えて、さらにその逆数の分布が使えて、その期待値が$\frac{1}{n-2}$であったので、上の確かめにおいて
$$
\begin{equation}
\frac{1}{A+a} = E(\frac{n-2}{|\mathbf{y}|^2})
\end{equation}
$$
と$n-2$を用いた。

今、$\mathbf{\theta} \sim N(M,A)$ であるとすると、$\frac{1}{A+a}$ の期待値は
$$
\begin{equation}
\frac{1}{A+a} = E(\frac{n-2}{|\mathbf{y-M}|^2})
\end{equation}
$$
である。これは変わらない。

しかしながら、Mが不明であるとすると、
$$
\begin{equation}
\hat{M}=\frac{\sum_{i=1}^n y_i}{n}
\end{equation}
$$
としたうえで、$\hat{M}$を$M$の代わりに使わなくてはならないが、自由度が1減っているから、
$$
\begin{equation}
\frac{1}{A+a} = E(\frac{n-2-1}{|\mathbf{y}-\hat{M}|^2})=E(\frac{n-3}{|\mathbf{y}-\hat{M}|^2})
\end{equation}
$$

のように$n-3$を使わなくてはならなくなる。

この$n-3$がよいことを確かめてみる。

```{r}
# 適当に真値と誤差の分散を決める
A <- 7
a <- 3
# シミュレーション回数
n.iter <- 100000
# Xはn-2を用いる、Yはn-3を用いる
# 真値が0を中心とした場合
X <- Y <- rep(0,n.iter)
# 真値が0でない値Mを中心とした場合で
# Mが既知である場合
X.M <- Y.M <- rep(0,n.iter)
# Mが未知である場合
X.M.2 <- Y.M.2 <- rep(0,n.iter)
M <- 20
for(i in 1:n.iter){
  n <- 100
  # 0,Mそれぞれを中心に真値を定める
	theta <- rnorm(n,0,sqrt(A))
	theta.M <- rnorm(n,M,sqrt(A))
  # 誤差はすべての場合で同一の値ベクトルを用いる
	z <- rnorm(n,0,sqrt(a))
  # 観測値
	y <- theta + z
	y.M <- theta.M + z
  # A/(A+a)の推定をそれぞれの方法で行う
	X[i] <- 1-(n-2)*a/sum(y^2)
	Y[i] <- 1-(n-3)*a/sum(y^2)
	X.M[i] <- 1-(n-2)*a/sum((y.M-M)^2)
	Y.M[i] <- 1-(n-3)*a/sum((y.M-M)^2)
	X.M.2[i] <- 1-(n-2)*a/sum((y.M-mean(y.M))^2)
	Y.M.2[i] <- 1-(n-3)*a/sum((y.M-mean(y.M))^2)
}
# 推定される値の答えと推定値を表示
A/(A+a)

mean(X)-A/(A+a)
mean(Y)-A/(A+a)

mean(X.M)-A/(A+a)
mean(Y.M)-A/(A+a)

mean(X.M.2)-A/(A+a)
mean(Y.M.2)-A/(A+a)

```
$n-2$を使ったmean(X)は$n-3$を使ったmean(Y)よりも良い推定値であり、真値の期待値がわかっているならば、$n-2$を使ったmean(X.M)が$n-3$を使ったmean(Y.M)よりも良いが、真値の期待値がわからず、その推定値を用いる場合には、$n-2$を使ったmean(X.M.2)は$n-3$を使ったmean(Y.M.2)よりも悪い様子がわかる。

従って、真値の平均値を観測値から推定するときには$n-3$を用いて係数を定めることがわかる。

次節でそのような関数を作っておこう。

## aがわかっていれば、もう推定できる。

ここまでで$M$の期待値と$\frac{1}{A+a}$の期待値が観測値$\mathbf{y}$のみから推定できた。

従って、$A$もしくは$a$を与えれば、$\frac{A}{A+a}$または$\frac{a}{A+a}$の期待値はわかり、真値の推定が可能となる。

誤差の分散aが既知であるとして、真値推定をしてみよう。

まずはそのための関数を以下のように作成する。

```{r}
js.3 <- function(y,a){
  n <- length(y)
  M <- mean(y)
  # n-2の代わりにn-3を使っている
  Inv.Aa <- (n-3)/sum((y-M)^2)
  k <- 1-a*Inv.Aa
  ret <- k*(y-M) + M
  return(ret)
}
```

真値も誤差も正規分布の場合を見てみる。

```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
m <- 0
A <- 5 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
Rs <- Rs.2 <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rnorm(n,M,sqrt(A))
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  theta.hat.2 <- js.3(y,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.2[i] <- sum((theta-theta.hat.2)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate, M,A known",ylab="ML estimate")
abline(0,1,col=2)
plot(Rs.2,Rs.ml,pch=20,asp=1,xlab="JS estimate, M,A unknown",ylab="ML estimate")
abline(0,1,col=2)
```

真のM,Aがわかっている方がよいが、わからなくても、良いことは両者に強い相関があることでわかる。

```{r}
plot(Rs,Rs.2,pch=20,asp=1,xlab="M,A known",ylab="M,A,unknown")
abline(0,1,col=2)
```

真値分布を指数分布にしてみる。
誤差は正規分布である。


```{r}
n.iter <- 1000
n <- 50
M <- 2 # 真値の平均
A <- M^2 # 真値の分散
m <- 0
a <- 7 # 誤差の分散を変更
Rs <- Rs.2 <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  theta <- rexp(n,1/M)
  z <- rnorm(n,m,sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  theta.hat.2 <- js.3(y,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.2[i] <- sum((theta-theta.hat.2)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate, M,A known",ylab="ML estimate")
abline(0,1,col=2)
plot(Rs.2,Rs.ml,pch=20,asp=1,xlab="JS estimate, M,A unknown",ylab="ML estimate")
abline(0,1,col=2)
plot(Rs,Rs.2,pch=20,asp=1,xlab="M,A known",ylab="M,A,unknown")
abline(0,1,col=2)
```

最後に誤差も正規分布ではなくしてみる。
```{r}
n.iter <- 1000
n <- 50
M <- -2 # 真値の平均を変更
A <- (M*2)^2/12 # 真値の分散を変更
a <- 7 # 誤差の分散を変更
Rs <- Rs.2 <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  # 平均Mの一様分布
  theta <- runif(n)*2*M
  # 分散がaの指数分布乱数
  z <- rexp(n,1/sqrt(a))
  y <- theta + z
  theta.hat <- js.1(y,M,A,a)
  theta.hat.2 <- js.3(y,a)
  Rs[i] <- sum((theta-theta.hat)^2)
  Rs.2[i] <- sum((theta-theta.hat.2)^2)
  Rs.ml[i] <- sum((theta-y)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate, M,A known",ylab="ML estimate")
abline(0,1,col=2)
plot(Rs.2,Rs.ml,pch=20,asp=1,xlab="JS estimate, M,A unknown",ylab="ML estimate")
abline(0,1,col=2)
plot(Rs,Rs.2,pch=20,asp=1,xlab="M,A known",ylab="M,A,unknown")
abline(0,1,col=2)
```

$\frac{1}{A+a}$の推定値を使うと、その真の値を使うよりはかなり悪いが、それでも、観測値そのものを使うよりはよいことがわかる。

## 誤差の分散 a を推定しながらやってみる

誤差の分散を与えられれば、複数観察による推定の改善が得られることがわかった。

では、どのようにして誤差の分散を与えればよいだろうか。

### 試験結果から、受験者の正答する確率(正答する力)を推定する

あるクラスにn人の学生がいるとする。
できる学生もできない学生もいる。

ある教科に関する多数の設問を正答する力というものを個々の学生が持っているものとする。
これを「正答確率の真値$\mathbf{\theta}=(\theta_1,...,\theta_n)$」とする。

今、N問の設問を答えさせたところ、N問中$\mathbf{y}=(y_1,...,y_n)$問ずつ正解した(これが観測値)。

さて個々の学生の$\theta_i$の推定値$\hat{\theta_i}$を出してみよう。

$\hat{\theta}_{ML} = \frac{y_i}{N}$という推定値があるだろう。

James-Stein流の推定値を出すには、誤差分散aが必要になる。

一般に、確率$\theta$で成功し、確率$1-\theta$で失敗するような事象を$N$回繰り返すと、成功回数は二項分布を取る。

二項分布の期待値は$N\theta$であり、その分散は$N\theta(1-\theta)$であることが知られている。

従って、$N$個の設問について成功確率$\theta_i$のもとで観察を行えば、成功確率の真値を中心とした$\frac{y_i}{N}$の分散は$\frac{N\theta(1-\theta)}{N^2}=\frac{\theta(1-\theta)}{N}$である。

確かめてみる。

```{r}
N <- 30
# 正答率として裾野の厚いベータ分布を使ってみる
theta <- rbeta(1,5,3)
n.iter <- 1000
x <- rep(0,n.iter)
for(i in 1:n.iter){
  tmp <- sample(0:1,N,replace=TRUE,prob=c(1-theta,theta))
  x[i] <- sum(tmp)/N
}
theta
mean(x)
theta*(1-theta)/N
var(x)
```

誤差分散は$\frac{\theta(1-\theta)}{N}$で良さそうである。

しかしながら残念なことに、$\theta$はn人についてばらばらであるし、そもそも$\theta$が推定の対象である。

今、欲しいのは、全員について等しく適用したい誤差分散である。

ここで$\theta$の期待値として、全員の成功確率の平均の期待値$\hat{\bar{\theta}}=\frac{\sum_{i=1}^n y_i}{n}$を使うことにしてはどうだろうか。

この値が、誤差分散としていかがなものか、シミュレーションしてみる。

```{r}
N <- 30
n <- 15
# 正答率として裾野の厚いベータ分布を使ってみる
thetas <- rbeta(n,5,3)
# 極端にできのよい学生と不出来な学生を作っておく
thetas[1] <- 0.98
thetas[2] <- 0.05
n.iter <- 1000
x <- matrix(0,n.iter,n)
y <- rep(0,n.iter)
for(i in 1:n.iter){
  for(j in 1:n){
    tmp <- sample(0:1,N,replace=TRUE,prob=c(1-thetas[j],thetas[j]))
    x[i,j] <- sum(tmp)/N
  }
  mean.theta <- mean(x[i,])
  y[i] <- mean.theta*(1-mean.theta)/N
}
z <- mean(thetas*(1-thetas)/N)
z
mean(y)
```

オーダーは外れていないが、とてもよい推定値になっているわけではない。

とはいえ、他にaとして使えるよい値もないので、これを使ってやってみよう。

```{r}
N <- 10
n <- 100
# 正答率として裾野の厚いベータ分布を使ってみる
thetas <- rbeta(n,5,3)
# 極端にできのよい学生と不出来な学生を作っておく
thetas[1] <- 0.98
thetas[2] <- 0.05
n.iter <- 1000
x <- matrix(0,n.iter,n)
y <- rep(0,n.iter)
# 真値からのずれ
Rs <- Rs.ml <- rep(0,n.iter)
for(i in 1:n.iter){
  for(j in 1:n){
    tmp <- sample(0:1,N,replace=TRUE,prob=c(1-thetas[j],thetas[j]))
    x[i,j] <- sum(tmp)/N
  }
  mean.theta <- mean(x[i,])
  y[i] <- mean.theta*(1-mean.theta)/N
  
  theta.ML <- x[i,]
  # 分散aとしてy[i]を指定する
  theta.JS <- js.3(x[i,],a=y[i])
  Rs[i] <- sum((thetas-theta.JS)^2)
  Rs.ml[i] <- sum((thetas-theta.ML)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1,xlab="JS estimate,M,A known",ylab="ML estimate")
abline(0,1,col=2)
```

James-Stein Estimatorの方がよいことがわかる。

繰り返し実験のうちの最後の推定結果を示せば以下の通り。
```{r}
plot(thetas,theta.ML,pch=20,xlim=c(0,1),ylim=c(0,1),xlab="True Values",ylab="Estimated Values")
points(thetas,theta.JS,pch=20,col=2)
```
James-Stein推定値(赤)は真値が低いときには少し高めに、真値が高いときには少し低めに出ている(誰もが全員の平均値の方に縮められている)。

試験の設問数が少ないと、実力はあるのにたまたま正答数が少ない学生も現れる。その学生の実力は、実際の正答割合より底上げしておく方が真の値に近い可能性があるし、たまたま正答数が多かった学生についても逆方向で同じことが言える。

ただし、そのようにすることによって、本当に実力がない学生、本当に実力がある学生の実力の推定値は真値よりも控えめになることになる。

そうなったとしても、全員について推定値と真値との乖離を乖離の２乗の和で評価する場合には、全員分の乖離の総量を小さくできる、とそういうわけである。

James-Stein 推定量を使うことに躊躇する理由の一つが、

> 本当は個人差があるのに、その個人差がないかのようにみんなを似たり寄ったりの値にすることは、推定としてどうなの？

という疑問だろう。「ずれの小さい推定値」がほしいという気持ちの他に、「個々のサンプルの違いをきちんと表した推定値」がほしいという気持ちが推定作業の中にはあることを表しているということである。

したがって、どんな推定値がどういう理由で欲しいかによって、James-Stein推定量を用いるのがよいかどうかは変わってきそうである。

しかしながら、「全て」についての真値からのずれを小さくしたいなら、James-Stein推定量は良いのである。

## 線形回帰に適用してみる

上記の「正答力」と同じような条件であるが、「普段から真面目である」かどうか(たとえば出席率など)という情報があったとするとどうなるだろうか。

「真面目」だからといって、「正答力」が高いとは限らないけれど、そこにはそれなりの関係があるのは事実である。

その「真面目度」と「真の正答確率」とが次のようになっているものとする。

```{r}
N <- 10
n <- 100
majimedo <- runif(n)
# 正答率として裾野の厚いベータ分布を使ってみる
thetas <- rep(0,n)
for(i in 1:n){
  thetas[i] <- rbeta(1,5*majimedo[i]*0.8+1,4*(0.9-majimedo[i]*0.8))
}
plot(majimedo,thetas,xlim=c(0,1),ylim=c(0,1),xlab="Majimedo",ylab="true probability")
```

ある程度の正の相関がある。

今、線形回帰で直線を引くことができる。

```{r}
lm.out <- lm(thetas ~ majimedo)
pred.lm <- predict(lm.out)
plot(majimedo,thetas,pch=20,xlim=c(0,1),ylim=c(0,1))
abline(lm.out[[1]][1],lm.out[[1]][2],col=2)
points(majimedo,pred.lm,pch=20,col=3)
```

ここで、ある真面目度の学生については、その回帰直線でその真面目度に関して指し示す値を「その学生にとっての平均値」とみなして、分散に応じて縮めてやることにすれば、真面目度情報がないときに、全員の平均に向かって縮めたことと対応した結果となるだろう。

実際、それをやってみる。

全員を同じ値に向かって縮めるのではなく、線形回帰の直線の相当値に縮める関数が以下である。

観測値から、全平均を推定しそこからの分散を計算したときには自由度1の分を考慮して$n-2$の代わりに$n-3$を用いた。

今回は、線形回帰直線について傾きと切片の２変数を推定しているので、その変数２個相当分を考慮して$n-4$を用いる。

```{r}
js.lm <- function(y,a,f){
  n <- length(y)
  lm.out <- lm(y ~ f)
  pred.lm <- predict(lm.out)
  # n-2の代わりにn-4を使っている
  Inv.Aa <- (n-4)/sum((y-pred.lm)^2)
  k <- 1-a*Inv.Aa
  ret <- k*(y-pred.lm) + pred.lm
  return(ret)
}
```

```{r}
N <- 10
n <- 100
majimedo <- runif(n)
# 正答率として裾野の厚いベータ分布を使ってみる
thetas <- rep(0,n)
for(i in 1:n){
  thetas[i] <- rbeta(1,5*majimedo[i]*0.8+1,4*(0.9-majimedo[i]*0.8))
}
n.iter <- 1000
x <- matrix(0,n.iter,n)
y <- rep(0,n.iter)
# 真値からのずれ
Rs <- Rs.ml <- Rs.lm <- rep(0,n.iter)
for(i in 1:n.iter){
  for(j in 1:n){
    tmp <- sample(0:1,N,replace=TRUE,prob=c(1-thetas[j],thetas[j]))
    x[i,j] <- sum(tmp)/N
  }
  mean.theta <- mean(x[i,])
  y[i] <- mean.theta*(1-mean.theta)/N
  
  theta.ML <- x[i,]
  theta.JS <- js.3(x[i,],y[i])
  theta.JS.lm <- js.lm(x[i,],y[i],majimedo)
  Rs[i] <- sum((thetas-theta.JS)^2)
  Rs.ml[i] <- sum((thetas-theta.ML)^2)
  Rs.lm[i] <- sum((thetas-theta.JS.lm)^2)
}
plot(Rs,Rs.ml,pch=20,asp=1)
abline(0,1,col=2)
plot(Rs.lm,Rs.ml,pch=20,asp=1)
abline(0,1,col=2)
```

全体の平均に向かって縮めても、線形回帰直線の値に向かって縮めても、観察値そのものを推定値とするよりも、全体のずれが小さいことがわかる。


全体の平均に向かうか、線形回帰に向かうかを比較すると、回帰直線に向かう方が当てはまりがよいことがわかる(縦軸が線形回帰直線に向かう方なので)。

```{r}
plot(Rs,Rs.lm,pch=20,asp=1)
abline(0,1,col=2)
```


繰り返し実験のうちの最後の推定結果を示せば以下の通りで、黒が真値、灰色が観察値(かつ、最尤推定値)、赤い直線が線形回帰直線。緑の点は全体平均へと縮めた結果で、青の点は、回帰直線へ向けて縮めた結果である。

全体平均への縮めた点(緑)では傾きが目立たなくなっているのに大して、回帰直線へ縮めた点(青)では、傾きを残したまま打点範囲が狭まっている様子が見て取れる。


```{r}
lm.out <- lm(theta.ML~majimedo)
pred.lm <- predict(lm.out)
plot(majimedo,thetas,pch=20,xlim=c(0,1),ylim=c(0,1),xlab="Majimedo",ylab="Thetas")
abline(lm.out[[1]][1],lm.out[[1]][2],col=2)
points(majimedo,theta.ML,pch=20,col="grey")
points(majimedo,theta.JS,pch=20,col=3)
points(majimedo,theta.JS.lm,pch=20,col=4)

```

縮みの様子を線で結べば以下のようになる。
値が大きい方に縮んだ場合に赤い線で、値が小さくなる方に縮む場合に黒い線で示している。


```{r}
lm.out <- lm(theta.ML~majimedo)
pred.lm <- predict(lm.out)
plot(majimedo,thetas,pch=20,xlim=c(0,1),ylim=c(0,1),xlab="Majimedo",ylab="Thetas")
abline(lm.out[[1]][1],lm.out[[1]][2],col=2)
points(majimedo,theta.ML,pch=20,col="grey")
points(majimedo,theta.JS,pch=20,col=3)
points(majimedo,theta.JS.lm,pch=20,col=4)
for(i in 1:n){
  col1 <- 1
  col2 <- 1
  if(theta.JS[i]-theta.ML[i]>0){
    col1 <- 2
  }
  if(theta.JS.lm[i]-theta.ML[i]>0){
    col2 <- 2
  }
  segments(majimedo[i],theta.ML[i],majimedo[i],theta.JS[i],col=col1)
  segments(majimedo[i],theta.ML[i],majimedo[i],theta.JS.lm[i],col=col2)
}
```
この縮み方の違いは次のようにも表現できる。

線形回帰直線へ縮めた場合には、直線の下の点はどんなに縮んでも直線の下の領域まで、直線の上の点は直線より上の領域の範囲内で縮んでいるのに対して、全体平均へと縮めた場合(緑)には、回帰直線を乗り越えて縮んだり、回帰直線から遠ざかるように縮む場合も現れる。

全体平均への縮みについてだけ打点し、線で結んでみる。
回帰直線より上であるのに、さらに上方に移動していたり、その逆に回帰直線より下方にあるのに、さらに下方に移動していることがあることが見て取れる。

これは、水平線で示した、観察値の全体平均に向かう縮みだからである。

```{r}
lm.out <- lm(theta.ML~majimedo)
pred.lm <- predict(lm.out)
plot(majimedo,thetas,pch=20,xlim=c(0,1),ylim=c(0,1),xlab="Majimedo",ylab="Thetas")
abline(lm.out[[1]][1],lm.out[[1]][2],col=2)
abline(h=mean(theta.ML),col=4)
points(majimedo,theta.ML,pch=20,col="grey")
points(majimedo,theta.JS,pch=20,col=3)
for(i in 1:n){
  col1 <- 1
  if(theta.JS[i]-theta.ML[i]>0){
    col1 <- 2
  }
  segments(majimedo[i],theta.ML[i],majimedo[i],theta.JS[i],col=col1)
}
```

それに対して、回帰直線への縮みでは、回帰直線に向かった縮みであることが見て取れる。

```{r}
lm.out <- lm(theta.ML~majimedo)
pred.lm <- predict(lm.out)
plot(majimedo,thetas,pch=20,xlim=c(0,1),ylim=c(0,1),xlab="Majimedo",ylab="Thetas")
abline(lm.out[[1]][1],lm.out[[1]][2],col=2)
points(majimedo,theta.ML,pch=20,col="grey")
points(majimedo,theta.JS.lm,pch=20,col=4)
for(i in 1:n){
  col2 <- 1
  if(theta.JS.lm[i]-theta.ML[i]>0){
    col2 <- 2
  }
  segments(majimedo[i],theta.ML[i],majimedo[i],theta.JS.lm[i],col=col2)
}
```

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------