当選確率nのくじをn回ひいた時に少なくとも1回は当る確率の計算について

先日このツイートを見て、念のためそれを確かめてみた。

ガチャの予算計算をするときに考える覚えることは3個だけあればいいと考えてる。
ひとつはさっき言った
・n分の1のガチャをn回引いて当たりが出る確率は63%
で、それに追加して
・3n回引くと95%で当たりが出る
・5n回引くと99%で当たりが出る
だけ覚えれば破産はしないと思う
— SCP-███ Lolida(ろりいだ) (@Lolida21) 2017年1月12日

n回中、少なくとも1回は当る確率

「n回の試行中、少なくとも1回は当たりが出る確率」は、「n回の試行中、1回も当たらない確率」を「1（＝100％）」から引けばいい。 *1

当る確率が「p」だとして、1回の試行で当たらない確率は「1- 1/p」。 n回連続で当たらない確率は、「（1-1/p）ⁿ」（累乗計算）となります。

試しにPythonで数式計算をしてみる

Pythonで数式計算をするライブラリは「sympy」

from sympy import *

sympyでは先ず使用する変数を定義します。ここでは1回の当選確率を「p」。求めたい値であるn回試行中1回は当る確率を「p1」。試行回数を「n」としました。

n,p,p1=symbols("n p p1")

先ほどの計算式を記述

p1 = 1 - (1-1/p)**n

また、今の設定では、当選確率の逆数の回数分だけくじを引くことになっているので「p == n」を代入し、それを新たに「p1_2」と定義します。

p1_2 = p1.subs(p,n)
print(p1_2)

すると、「1/n」の確率で当るくじを「n」回ひいた時に、少なくとも1回は当るが出る確率は

「-(1 - 1/n)**n + 1」

とわかりました。先ほどの手計算の数式と同じです。

このp1_2に具体的に「n==2」なり「n==100」なりを代入した結果を出力すれば、求めたい数値が出てきます。「evalf()」で分数が小数点表示になります

print(p1_2.subs(n,2).evalf())
print(p1_2.subs(n,6).evalf())
print(p1_2.subs(n,100).evalf())

当選確率1/2で2回だと「0.750000000000000」 1/6の確率で6回だと「0.665102023319616」 1/100の確率で100回だと「0.633967658726771」

値を求めるだけならRでできる

数式を求めることに意義が薄いと思えば最初からRを使えば手っ取り早いです。

当たり／外れ　の2種類しかない試行は二項分布になります。確率を求めるのは、「d～」で、二項分布は「binom」。引数は、「当たりの確率」、「試行回数」、「当る確率」の3つで、先ほど同様、1回も当たらない（＝＝０）確率を1から引けば良いです。

1-dbinom(0,2,1/2)
1-dbinom(0,4,1/4)
1-dbinom(0,6,1/6)
1-dbinom(0,100,1/100)

> 1-dbinom(0,2,1/2)
[1] 0.75
> 1-dbinom(0,4,1/4)
[1] 0.6835938
> 1-dbinom(0,6,1/6)
[1] 0.665102
> 1-dbinom(0,100,1/100)
[1] 0.6339677

それはそうと、さっきと数字が違うじゃないですか？

さっきのツイートの数字は、63％と断言されていましたが、実際には違う、、、でも1/100で100回の時はとても近い

Rで繰り返し処理

というわけで、n==2　から　n==100　までの夫々の確率をデータフレームに格納して、プロットしてみます。

繰り返し処理でデータフレームに追記していく場合、rbindを使います。

df <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n,1/n))
  df <- rbind(df,dat)
}
colnames(df) <- c("prob")
df$x <- rep(2:100)

結果を出力すると、試行回数が上がるに連れて一定の数字に修練しているようです。

> head(df)
       prob x
1 0.7500000 2
2 0.7037037 3
3 0.6835938 4
4 0.6723200 5
5 0.6651020 6
6 0.6600833 7
> tail(df)
        prob   x
94 0.6340653  95
95 0.6340450  96
96 0.6340250  97
97 0.6340055  98
98 0.6339864  99
99 0.6339677 100

グラフにしてみる

ggpotでグラフ化

p <- ggplot(data=df,aes(x=x,y=prob))
p <- p + geom_bar(stat="identity")
p <- p + ggtitle("当たりの確率の逆数と同じ回数だけひいたら１回は辺りが出る確率")
p <- p + geom_text(aes(x=100,
                       y=df[df$x==100,]$prob,
                       label=round(df[df$x==100,]$prob*100,2)),
                   vjust=-1)
plot(p)
ggsave(plot=p,file="20170115.png",h=4,w=6,scale=2)

ほぼほぼ63.4%に収斂しています。

f:id:yyhhyy:20170115172459p:plain

それでは、「3n回」や「5n回」は？

繰り返し計算中の確率計算部分で試行回数を「n3」や「n5」にするだけです

dat <- data.frame(1- dbinom(0,n*3,1/n))

3倍の回数なら

f:id:yyhhyy:20170115172557p:plain

5倍の回数なら

f:id:yyhhyy:20170115172609p:plain

実際のところ、確率1％であっても100回も引く人はいないと思いますが、、、

これだけ回数を増やせばポアソン分布である

ガチャの確率が一般的にどれくらいなのか知りませんが、二項分布は確率が低く、試行回数が増えればポアソン分布で近似できる、ということがわかっています。

二項分布 - Wikipedia

dskjal.com

qiita.com

例えば、当選確率pが「1/100」で試行回数nが「100回」なのであれば、

lambda = n * p は、「100*(1/100)」

同様に「300回」なら「100(1/300)」、「500回」なら「500(1/100)」を引数として

1 - dpois(0,(100*(1/100)))
1 - dpois(0,(300*(1/100)))
1 - dpois(0,(500*(1/100)))

とすれば、先ほど収束していった値が出てきます。

> 1 - dpois(0,(100*(1/100)))
[1] 0.6321206
> 1 - dpois(0,(300*(1/100)))
[1] 0.9502129
> 1 - dpois(0,(500*(1/100)))
[1] 0.9932621

コード全体

from sympy import *
n,p,p1=symbols("n p p1")
p1 = 1 - (1-1/p)**n
p1_2 = p1.subs(p,n)
print(p1_2)
print(p1_ans.evalf())
print(p1_2.subs(n,6).evalf())
print(p1_2.subs(n,100).evalf())

library("dplyr", lib.loc="C:/hogehoge/r/R-3.3.2/library")
library("reshape2", lib.loc="C:/hogehoge/r/R-3.3.2/library")
library("ggplot2", lib.loc="C:/hogehoge/r/R-3.3.2/library")

1-dbinom(0,2,1/2)
1-dbinom(0,4,1/4)
1-dbinom(0,6,1/6)
1-dbinom(0,100,1/100)

df <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n,1/n))
  df <- rbind(df,dat)
}
colnames(df) <- c("prob")
df$x <- rep(2:100)
head(df)
tail(df)

p <- ggplot(data=df,aes(x=x,y=prob))
p <- p + geom_bar(stat="identity")
p <- p + ggtitle("当たりの確率の逆数と同じ回数だけひいたら１回は辺りが出る確率")
p <- p + geom_text(aes(x=100,
                       y=df[df$x==100,]$prob,
                       label=round(df[df$x==100,]$prob*100,2)),
                   vjust=-1)
plot(p)
ggsave(plot=p,file="20170115.png",h=4,w=6,scale=1)

df_3 <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n*3,1/n))
  df_3 <- rbind(df_3,dat)
}
colnames(df_3) <- c("prob")
df_3$x <- rep(2:100)
print(df_3)

p3 <- ggplot(data=df_3,aes(x=x,y=prob))
p3 <- p3 + geom_bar(stat="identity")
p3 <- p3 + ggtitle("当たりの確率の逆数との３倍の回数をひいたら１回は辺りが出る確率")
p3 <- p3 + geom_text(aes(x=100,
                       y=df_3[df_3$x==100,]$prob,
                       label=round(df_3[df_3$x==100,]$prob*100,2)),
                   vjust=-1)
plot(p3)
ggsave(plot=p3,file="20170115_3.png",h=4,w=6,scale=1)


df_5 <- data.frame()
for (n in 2:100){
  dat <- data.frame(1- dbinom(0,n*5,1/n))
  df_5 <- rbind(df_5,dat)
}
colnames(df_5) <- c("prob")
df_5$x <- rep(2:100)
print(df_5)

p5 <- ggplot(data=df_5,aes(x=x,y=prob))
p5 <- p5 + geom_bar(stat="identity")
p5 <- p5 + ggtitle("当たりの確率の逆数との３倍の回数をひいたら１回は辺りが出る確率")
p5 <- p5 + geom_text(aes(x=100,
                         y=df_5[df_5$x==100,]$prob,
                         label=round(df_5[df_5$x==100,]$prob*100,2)),
                     vjust=-0.5)
plot(p5)
ggsave(plot=p5,file="20170115_5.png",h=4,w=6,scale=1)

1 - dpois(0,(100*(1/100)))
1 - dpois(0,(300*(1/100)))
1 - dpois(0,(500*(1/100)))