合格率60% の認定試験でほとんどの国民は一般人になれるか

MikuHatsune2017-06-19

虚構新聞のネタなので、一般人を試験により認定するという話は存在しません。
 
こんな記事を見かけた。
一般人認定試験、来年度実施を検討 「共謀罪」成立受け
共謀罪について、一般人であるかどうかの認定試験を行って、一般人か否かが決まる。ここで、「合格率は60%程度。繰り返し受検すればほとんどの人が一般人になれる」と言っている。
 
ここで、合格率が60% の試験を何回受けたらほとんどの人が合格して一般人になれるか考える。
10万人の受検者がいて、1回の試験の合格率が60% のとき、このシミュレーションでは13回目で全員が合格した。

n <- 100000
p <- 0.6
x <- 0:15
res <- c(n, rep(0, length(x)-1))
for(i in 1:(length(res)-1)){
  m <- rbinom(res[i], 1, p)
  tab <- table(factor(m, 0:1))
  res[i+1] <- tab[1] # 合格していない人
}

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, las=1)
plot(head(x+1, -1), tail(res/n, -1), xlab="初めて合格したときに受検していた回数", ylab="合格できていない受験者の割合")
 [1] 100000  40084  16161   6473   2627   1057    424    166     56     25      5      1      1      0
 [15]      0      0

 
各施行がベルヌーイ施行のとき、1回の試験での0/1 binary アウトカムは二項分布から得られた。ここで、二項分布の拡張で負の二項分布というものがある。これは、r 回の成功を得るまでに必要な施行回数、とか、r 回の成功を刷る前に失敗した試行回数、とか言われる。
ここで、1回目にいきなり成功する確率と、1回目は失敗して2回目に成功する確率の和を考える。成功確率が0.6 のとき、r+(1-r)r で与えられる。
欲しい成功回数が1回のときは、負の二項分布は幾何分布 *geom と同じである。

p <- 0.6
p + (1-p)*p
[1] 0.84

ここで、今回の状況では、1回の成功(合格)までに何回受検するかだから、dnbinom では

dnbinom(x, 1, p)
 [1] 6.000000e-01 2.400000e-01 9.600000e-02 3.840000e-02 1.536000e-02 6.144000e-03 2.457600e-03
 [8] 9.830400e-04 3.932160e-04 1.572864e-04 6.291456e-05 2.516582e-05 1.006633e-05 4.026532e-06
[15] 1.610613e-06 6.442451e-07

pnbinom では

pnbinom(x, 1, p)
 [1] 0.6000000 0.8400000 0.9360000 0.9744000 0.9897600 0.9959040 0.9983616 0.9993446 0.9997379 0.9998951
[11] 0.9999581 0.9999832 0.9999933 0.9999973 0.9999989 0.9999996

となる。ここで、x は0 から始めているので、1回の成功までに失敗する回数が返ってくる。
 
理論的に求めると、10回受検しても10^{-4}、つまり1万人に1人は合格できないし、10万人に1人にするなら13回必要になる。

h <- pnbinom(x, 1, p)
e10 <- 10^(-(2:5))

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, las=1)
plot(x+1, log10(1-h), xlab="初めて合格したときに受検していた回数", ylab="合格できていない受験者の割合 [log10]")
abline(h=log10(e10), lty=3)


 
シミュレーションと重ねてみると、受検者が減った10回目くらいまでの試験では理論値とよく一致する。

par(mar=c(4.5, 5, 2, 2), cex.lab=1.5, las=1)
plot(x+1, log10(1-h), xlab="初めて合格したときに受検していた回数", ylab="合格できていない受験者の割合 [log10]")
points(head(x+1, -1), tail(log10(res/n), -1), pch=15)
abline(h=log10(e10), lty=3)