r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

ババ抜きで13枚以上のスタートになる確率

ババ抜きで13枚スタートになる確率

異動が決まり、週末は長距離移動かあるいはフルに仕事。
サンデープログラマーなので最近記事がかけていない。
秋からは精力的に更新したいものだ。

某アニメでババ抜きのシーンが有った

4人でゲームを行い、なんと13枚スタートであった。
これは、相当にレアなのではないか?
Rで計算&シミュレートしてみる

ある個人が13枚配られた時に13枚スタートになる確率を計算する。
勝率をコントロールするため、自分で配って有利な枚数のものを取る
輩がいなければ、この結果を3/4倍すれば求める確率になる。

13枚以上スタートになる確率だと14枚スタートになる確率を
足してやる必要がある。

まずは13枚配られた時に13枚スタートになる確率

#小数で
library(Rmpfr)
fact1 <- function(x){
      one <- mpfr(1,1000)
      fact.x <- Reduce("*", c(one, 1:x))
}

four <- mpfr(4,1000)
a13 <- (four ^ 13 + 13 * four ^ 12) * fact1(13) * fact1(40) / fact1(53)
a13

#分数で

library(gmp)
fact2 <- function(x){
      one <- as.bigz(1)
      fact.x <- Reduce("*", c(one, 1:x))
}

four2 <- as.bigz(4)
a13.2 <- (four2 ^ 13 + 13 * four2 ^ 12) * fact2(13) * fact2(40) / fact2(53)
a13.2

小数
0.0003389767722882 以下略

分数
8388608/24746851955

2950.05464017391204833984375回に1回起こる。
約3000回に一回なので相当レアだ。

シミュレーションで確認する

tra <- c(rep(1:13, 4), 14)

sim = 1000000
res1 <- numeric(sim)
for(i in 1:sim){
      temp <- sample(tra, 13, replace = FALSE)
      temp2 <- table(temp)
      res1[i] <- length(temp2[temp2 == 1 | temp2 == 3])
}
sum(res1 > 12) / sim
table(res1)
hist(res1, freq = FALSE, xlim = c(0,14))

> sum(res1 > 12) / sim
[1] 0.000332

0.0003389767722882 と極めて近いのでOKだろう。

> binom.test(sum(res1 > 12) , sim)

95 percent confidence interval:
0.0002972492 0.0003696965

0.0003389767722882を含んでいる。

次に14枚スタートになる確率

a14 <- four ^ 13 * fact1(13) * fact1(40) / fact1(53)
a14
1 / a14

> ans14
[1] 7.97592405以下略e-5

> 1 / ans14
[1] 12537.73222073以下略

13枚以上スタートになる確率は

ans <- a13 * 3 / 4 + a14 * 1 / 4
ans
1 /ans

> ans <- a13 * 3 / 4 + a14 * 1 / 4
> ans
[1] 0.00027417238935以下略

> 1 /ans
[1] 3647.34028以下略

ある個人が13枚以上スタートになるのは、およそ3650回に一回のレアな現象であった。


4人の誰か1人以上に起こる、ということだと

ans2 <- 1 - (1 - a13)^3 * (1 - a14)
ans2
1/ans2

> ans2
[1] 0.00109626以下略

> 1/ans2

[1] 912.189202以下略

900回に一回程度であった。