読者です 読者をやめる 読者になる 読者になる

広告/統計/アニメ/映画 等に関するブログ

広告/統計/アニメ/映画 等に関するブログ

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

R python 統計 雑感

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

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%に収斂しています。

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("当たりの確率の逆数と同じ回数だけひいたら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)))

*1:数学が面白いのはこの手の逆転の発想が普通に出てくることですよね