先日このツイートを見て、念のためそれを確かめてみた。
n回中、少なくとも1回は当る確率
「n回の試行中、少なくとも1回は当たりが出る確率」は、「n回の試行中、1回も当たらない確率」を「1(=100%)」から引けばいい。
*1
当る確率が「p」だとして、1回の試行で当たらない確率は「1- 1/p」。
n回連続で当たらない確率は、「(1-1/p)n」(累乗計算)となります。
試しに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回も当たらない(==0)確率を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("当たりの確率の逆数と同じ回数だけひいたら1回は辺りが出る確率")
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%に収斂しています。
それでは、「3n回」や「5n回」は?
繰り返し計算中の確率計算部分で試行回数を「n3」や「n5」にするだけです
dat <- data.frame(1- dbinom(0,n*3,1/n))
3倍の回数なら
5倍の回数なら
実際のところ、確率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("当たりの確率の逆数と同じ回数だけひいたら1回は辺りが出る確率")
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("当たりの確率の逆数との3倍の回数をひいたら1回は辺りが出る確率")
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("当たりの確率の逆数との3倍の回数をひいたら1回は辺りが出る確率")
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)))