『確率４万分の１、県の入札くじに６回連続当選「奇跡的」』なのかどうかを調べよう

www.asahi.com

動機

面白いニュースを見かけたので計算してみましたなポストです。記事のこの部分に注目しました。

６回のくじにはそれぞれ制限価格で入札した３～８の鑑定業者が参加した。６回連続で当たる確率を計算すると約４万分の１となる。

奇跡ってどのくらいの割合で起きるのかな？もうすぐクリスマスだしってことでRでやってみます。

まず引用箇所を表にすると、

i回目	鑑定業者数 X_i
1	X_1
2	X_2
3	X_3
4	X_4
5	X_5
6	X_6

ここでX_i ∈ [3,8]の整数で、「６回連続で当たる確率を計算すると約４万分の１」を満たすX_iであって欲しいってことでX_iを全部掛け算したら4万になれば良いでしょう。

とりあえずそのような鑑定業者数の組み合わせを力技で探すと、

kgyousya <- c(3,4,5,6,7,8)

kaisuu <- data.frame(
  x1=kgyousya
  ,x2=kgyousya
  ,x3=kgyousya
  ,x4=kgyousya
  ,x5=kgyousya
  ,x6=kgyousya
)

kaisuu.d <- expand.grid(kaisuu)
kaisuu.m <- as.matrix(kaisuu.d)

kaisuu.prod <- apply(kaisuu.m,MARGIN = 1,prod)

table( kaisuu.prod[(kaisuu.prod > 39000) & (kaisuu.prod < 41000)] )

#> table( kaisuu.prod[(kaisuu.prod > 39000) & (kaisuu.prod < 41000)] )
#
#39200 40000 40320 40960 
#180    15   720    60

915ケースありました。

次に、掛け算して39,200〜40,960となるような鑑定業者数の組み合わせ：915ケースで、 各ケースを10000回シミュレーションする。 そして各ケースで確率4万分の1以下が10000回中何回起きるかを調べる。(←奇跡の割合） 最後に奇跡の割合の分布を観察しましょう。

# 候補のケース数: 915
#> nrow(kaisuu.sub)
#[1] 915
mycase <- nrow(kaisuu.sub)
results <- numeric(mycase)

# シミュレーション回数
ntimes <- 1:10000
result <- numeric( length(ntimes) )

for(j in 1:mycase){
  for(i in ntimes){
    np <- sample(c(0,1),size=6,replace = TRUE, prob=c(.5,.5))
    result[i] <- prod( abs( np  -  (1/kaisuu.sub[j,]) ) )
  }
  results[j] <- sum( result < 1/40000 )/length(result)
}

結果

度数分布表

> table(results)
results
    0 0.005 0.006 0.007 0.008 0.009  0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 
  180     1     3     7     8    17    21    56    63    66    64    54    74    74    49 
0.019  0.02 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 0.029 
   46    45    30    18    19     7     4     4     3     1     1 

要約値

> summary(results)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.01000 0.01400 0.01223 0.01700 0.02900

どんなケースでも平均すれば1%ぐらいは起きるのでは？という見方で良いのかなぁ・・・。

プロット

感想

こういうの式だけ立てて計算できるようになりたいですなぁ。