2014-03-08

実力を知るには仲間が必要

James-Stein

一昨日(こちら)、たくさんの受験者の結果を全部使うと、使わない場合より、個々人の真の正答力を、より正確に推定できるという話を書いた
より正確に推定する方法をJames-Stein Estimatorと言った
「ふーん」と思いやすい例はないかと考えた
こんな風に考える
今、ある科目について、設問に正解する確率というものを個人が持っていて、それを「真の正解力」と呼ぶことにする。
『私』がN問の試験を受けて、 $x_0$ 問正解したとする
『私』の真の正解力は $\frac{x_0}{N}$ と推定するのがよいだろう
この試験を一人きりで受けていたら、これで話は終わりである。
もし、一緒に受けていた友人がいて、その友人が「僕は $x_1$ 問正解していたよ」と教えてくれたとしよう
このとき、『私』の「真の正解力」はあいかわらず $\frac{x_0}{N}$ と推定するのがよい
２人目の友人が「 $x_2$ 問正解したよ」と教えてくれると $\frac{x_0}{N}$ が『私』の真の正解力の推定値としておくのがよいかどうかに影響をうけはじめる。『私』の「真の正解力」は $\frac{x_0}{N}$ と考えるよりも、自分と友人とのみんなの平均点を考慮する方が「真の正解力」に近いことが多いことが知られている、というのがJames-Stein推定です
友達の多い『私』なので、３人目の友人が「 $x_3$ だった」と教えてくれたとする。さらに影響をうける。
「『私』は友人と仲良しかもしれないが、正解力に関して関係ない」と思っても、完全に外れている場合(実は、自分はエイリアンだった、というくらいのはずれ具合)以外は、友人の正解数を参考にした方がよい、とこういうことである
その様子を友人数を1人ずつ増やしつつ、「友人を考慮した推定値」が「真の正解力」に近づく様子をグラフにしてみる
もちろん、１回の試験なので、どんなに友人数が増えても「真の正解力」に収束するわけではなく、１回だけの試験での『まぐれ』を発揮すれば、「真の正解力」をある程度隠ぺいできる
ただし注意が必要
「友人」の結果を参考することで推定精度が上がる、というのは、すべての人について見渡した場合であって、特定の個人についての推定結果が必ずよくなるというわけではない
詳細は一昨日の記事(こちら)で

# 二項分布の場合のJames-Stein推定
js.binom <- function(y,N){
  n <- length(y)
  x <- y/N
  M <- mean(x)
  # n-2の代わりにn-3を使っている
  a <- M*(1-M)/N
  Inv.Aa <- (n-3)/sum((x-M)^2)
  k <- 1-a*Inv.Aa
  ret <- k*(x-M) + M
  return(ret)
}
# 最大友人数
max.n <- 1000
# 友人の真の正答力
thetas <-rbeta(max.n,1,2)
# 『私』の真の正答力
my.theta <- 0.7

# 試験の設問数
Ns <- 1:100

# 最尤推定とＪＳ推定の結果を格納する行列
est.js <- est.ml <- matrix(0,length(Ns),max.n+1)

# 試験の実施結果の格納行列
test.result <- matrix(0,max.n+1,max(Ns))
# 『私』の試験結果
test.result[1,] <-cumsum(sample(0:1,max(Ns),replace=TRUE,prob=c(1-my.theta,my.theta)))
# 友人の試験結果
for(i in 1:max.n){
	test.result[i+1,] <-cumsum(sample(0:1,max(Ns),replace=TRUE,prob=c(1-thetas[i],thetas[i])))
}

# j問の試験だったとして、友人i-1人だったときの推定値を計算して格納

EST.ML <- EST.JS <- matrix(0,max.n+1,max(Ns))

for(i in 1:(max.n+1)){
	for(j in 1:max(Ns)){
		EST.ML[i,j] <- test.result[1,j]/j
		if(i<3){
			EST.JS[i,j] <- EST.ML[i,j]
		}else{
			tmp <- js.binom(test.result[1:i,j],j)
			EST.JS[i,j] <- tmp[1]
		}
		
	}
}

# 設問数がcheck.Nだった場合を表示
check.N <- 40
#par(mfcol=c(1,2))
# 友人は0人からn.friend人までを描図
n.friend <- 100
matplot(0:n.friend,cbind(EST.ML[1:(n.friend+1),check.N],EST.JS[1:(n.friend+1),check.N]),type="l",ylim = range(c(my.theta,EST.ML[,check.N],EST.JS[,check.N])),xlab="友人数",ylab="得点力・推定値")
abline(h=my.theta,col=3)
legend(50, 0.74, c("得点率そのまま","JS推定値","真の得点力"), col = 1:3, lty = c(1,2,1))
head(EST.JS[,check.N])