STAP細胞が見つかるまでひたすら適当な細胞を生成し続けるプログラムを見てお勉強

STAP細胞が見つかるまでひたすら適当な細胞を生成し続けるプログラム - 驚異のアニヲタ社会復帰への道

理論的な解が存在するかどうかをまず確認 - 裏 RjpWiki

いつもお勉強させて頂いているblogで面白そうなことをやっていた。

コードを参考にしながら、自分でもやってみる。

４つ組でsampleしたり、１つ組でsampleしたり。

確率密度関数を書いてみたりしてみた。

p = 1/26^4
k = 1
r = rnbinom(n = 1e+07, size = k, prob = p)

qnbinom(c(0.01, 0.05, 0.5, 0.95, 0.99, 0.9999, 0.999999), size = k, prob = p)

## [1]    4592   23439  316754 1368976 2104449 4208915 6313379

curve(dnbinom(x, size = k, prob = p), from = 0, to = 1e+07)


f0 <- function(){
stap <- "STAP"
tmp <- ""
n <- 0
while(tmp != stap){
 tmp <- paste(sample(LETTERS, size=4, replace=TRUE), collapse="")
 n <- n + 1
}
n
}

f1 <- function(){
jikken <- function(){sample(LETTERS, size=1, replace=TRUE)}

temp <- character(4)
n <- 0

temp[1] <- jikken()
temp[2] <- jikken()
temp[3] <- jikken()
temp[4] <- jikken()

while(temp[1] != "S" | temp[2] != "T" | temp[3] != "A" | temp[4] != "P") {
 temp[1:3] <- temp[2:4]
 temp[4] <- jikken()
 n <- n + 1
}
n
}

f2 <- function(){
jikken2 <- function(){sample(LETTERS, size=4, replace=TRUE)}
n <- 0
temp <- character(4)
while(temp[1] != "S" | temp[2] != "T" | temp[3] != "A" | temp[4] != "P") {
 temp <- jikken2()
 n <- n + 1
}
n
}

f0.t <- function(){
t1 <- proc.time()["elapsed"]
n0 <- f0()
t2 <- proc.time()["elapsed"]
as.numeric((t2 - t1)/n0 * 100000)
}

f1.t <- function(){
t1 <- proc.time()["elapsed"]
n1 <- f1()
t2 <- proc.time()["elapsed"]
as.numeric((t2 - t1)/n1 * 100000)
}

f2.t <- function(){
t1 <- proc.time()["elapsed"]
n2 <- f2()
t2 <- proc.time()["elapsed"]
as.numeric((t2 - t1)/n2 * 100000)
}

f0.time <- replicate(20, f0.t())
f1.time <- replicate(20, f1.t())
f2.time <- replicate(20, f2.t())

boxplot(f0.time, f1.time, f2.time)
t.test(f0.time, f1.time) #各p値は脳内で３倍する
t.test(f0.time, f2.time)
t.test(f1.time, f2.time)
summary(f0.time)
summary(f1.time)
summary(f2.time)

 > summary(f0.time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.057 1.300 1.313 1.305 1.335 1.521
> summary(f1.time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.691 1.712 1.716 1.729 1.725 1.899
> summary(f2.time)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.059 1.065 1.074 1.074 1.079 1.101

単位sampleあたりは、f2が最速のようだ。