2013-08-11

Chapter 1 Empirical Bayes and the James-Stein Estimator ぱらぱらめくる『Large-Scale Simultaneous Inference (講義＠Stanford)』

考え方は1950年代に遡るが、同時並行でたくさんのデータ解析を行うようになって実問題に適用されるにいたった。その最初の適用対象がマイクロアレイ・体系的発現解析
(ベイズ)推定なのか、検定なのかの区別なども問題になる
1.1 ベイズ規則と多変量正規分布推定
- N個が正規分布からのサンプルであるときに、それを同時に観察したときに、それぞれの観測値から、それぞれの元の値を推定しましょう、という話
  - - ここでは長さNのベクトル $\mu,\mathbf{0}$ , $N \times N$ 行列 $I$ を使って表している
- 平均0,分散Aの正規分布に従っている $\mu$ からサンプリングして、平均0分散1の正規乱数誤差で観測してzを得るとき、zの値から元の $\mu$ を推定すると、平均 $Bz$ ,分散 $B$ ( $B=\frac{A}{A+1}$ )の正規分布になる

# 試行回数(同時観測サンプル数)
N <- 10000
# 適当にAを取る
A <- runif(1)*10
# Bが求まる
B <- A/(A+1)
# サンプリングする
mu <- rnorm(N,0,sqrt(A))
# muは観測されないでこれに誤差を入れたzが観察される
z <- mu + rnorm(N,0,1)
#plot(mu,z)
# 観察されたzの値ごとに平均Bz,分散Bの正規分布から乱数を発生させてみる
# 標準正規乱数を取り
w <- rnorm(N,0,1)
# 分散をBにして
w <- w*sqrt(B)
# 平均をBzにする
w <- w + B*z
# 元のmuと、推定したwとを横軸に、観察値zを縦軸に重ねてプロット
plot(mu,z,xlim=range(c(mu,w)),ylim=range(z),cex=0.1,pch=20)
points(w,z,col=2,cex=0.1,pch=20)

- 最尤推定値は $z$ そのものであるのに対し、ベイズ推定値は $Bz$
- どちらの推定がよいかは、推定値と真値の差の二乗を全サンプルの和で評価してもよいだろう(リスク関数をそのように定める、ということ)
  - 最尤推定値の場合は、個々のサンプルで観測誤差を分散１の正規分布としているから、リスク関数の値(の期待値)は $N$
  - ベイズ推定の方は $N\frac{A}{A+1}$ になる
  - ベイズ推定の方が小さくなる…ただし、 $A$ がいくつかわからなければ $B$ もわからないわけで、 $A$ を事前情報として与える、という条件付きであるが
1.2 Empirical Bayes Estimation
- $A$ を知りたい
- たくさんのデータがあれば、推定はできる
- 実際、 $\hat{\mu} = Bz = \frac{A}{A+1}z=(1-\frac{1}{A+1})z$ なので、 $\frac{1}{A+1}$ の推定値がほしい
- 観察値のばらつきは、真値のばらつきと観察誤差のばらつきの和になっているから、観察値のばらつき具合から真値のばらつきを推定してやる
  - 正規分布を仮定しているから、観察値の普遍分散を出して…とやってもよさそうなところだが、テキストでは、二乗和を出して、それが自由度Nのカイ二乗分布様に従うことを使って推定している
- 実際、の推定値はとなる
  - 念のため確かめておこう

A <- runif(1)*10
N <- 100

n.trial <- 10000
SS <- rep(0,n.trial)
for(i in 1:n.trial){

mu <- rnorm(N,0,sqrt(A))
z <- mu + rnorm(N,0,1)

S <- sum(z^2)
SS[i] <- S
A.mean[i] <- S/N -1
A.median[i] <- S/(N-2)-1
}

A+1
mean(SS/N)

1/(A+1)
mean((N-2)/SS)

A+1
# Nで割った時が「よい」推定値
mean(SS/N)
mean(SS/(N-1))
mean(SS/(N-2))

# N-2を割った時が「よい」推定値
1/(A+1)
mean(N/SS)
mean((N-1)/SS)
mean((N-2)/SS)

> A+1
[1] 7.734166
> # Nで割った時が「よい」推定値
> mean(SS/N)
[1] 7.731988
> mean(SS/(N-1))
[1] 7.810089
> mean(SS/(N-2))
[1] 7.889784
> 
> # N-2を割った時が「よい」推定値
> 1/(A+1)
[1] 0.1292964
> mean(N/SS)
[1] 0.1319799
> mean((N-1)/SS)
[1] 0.1306601
> mean((N-2)/SS)
[1] 0.1293403

- これを使ったがJames-Stein 推定量
  - JS推定量を使って、真値からのずれを評価した値(リスク関数の値)を計算すれば、「真のA」が与えられた場合よりは大きくなる(けれど、それは仕方ない)
- 共通の分布からのサンプルについて、たくさん( $N \ge 3$ )のテストがあるときは個々の最尤推定量よりJS推定量を採用する方が、「全体のずれ」を小さくできる
- ここまでの議論では $\mu$ が0としてきたけれど、そうでない場合にもあてはまることを示すのは簡単

my.JS.estimate <- function(x){
	m <- mean(x) # ここに推定が入っているのでlength(x)-2でなく -3
	＃A1 <- (length(x)-2)/sum((x-m)^2)
	A1 <- (length(x)-3)/sum((x-m)^2)
	y <- (1-A1)*(x-m) +m
	return(y)
}

1.3 Estimating the Individual Components
- こんなにすごいのに、それほど使用されてこなかったのは、当てはまりの良しあしを「全て」の総合計で判断する方法であるので、個々のずれが大きくなる場合が出たことに(も)よる。外れ値がそれに影響する。
- 結構適当なデータに対してJS推定値を出してみよう

# 正規分布と指数分布の合わさったものをつくってみる
N <- 100
N2 <- 50
A <- 4
A2 <- 6
mu <- rnorm(N,3,sqrt(A))
z <- mu + rnorm(N,0,1)
mu2 <- rexp(N,1)+3
z2 <- mu2 + rnorm(N,0,2)
mus <- c(mu,mu2)
zs <- c(z,z2)
# JS推定
JSs <- my.JS.estimate(zs)
# リスク関数の値を見よう
sum((zs-mus)^2)
sum((JSs-mus)^2)
plot(mus,zs)
points(mus,JSs,col=2)

1.4 Learning from the Experience of Others
- ある特定のサンプルの推定をするに際して、他のサンプルのデータを使いながら(それに基づいた推定をする)、という意味でベイズ
- 他のサンプルのすべてを利用するか一部を利用するか、一部にするなら、どうやって選ぶか、ということを問題にすることもできる
1.5 Empirical Bayes Confidence Intervals
- Bの推定値に区間推定を持ちこむ、という手がある