r-statistics-fanの日記

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

レアなタクシーに遭遇する確率

レアなタクシーに遭遇する確率

京都で学会があり、レアなタクシーに遭遇する機会があった。

ヤサカタクシーは通常三つ葉なのだが、四つ葉と双葉
が存在する。それぞれ1200台中4台と2台しか無いらしい。

今回同じ日に四つ葉と双葉を見た。なんだか良い気分。

多分20-30台のタクシーを見たので、その場合に
四つ葉と双葉を両方見るのがどの位レアな事象なのか
気になった。

Rで計算してみる

四つ葉1台と双葉1台自身とそれ以上にレアな組合せ
すべての確率の合計を計算する。
同じタクシーに複数回遭遇することもあるという前提
で各遭遇は独立とする。

library(gmp)

f0 <- function(m=30, x=1, y=1){
n = 1200 #n 総台数
#m = m #m 出会った台数
a = 2 #a ふたば
b = 4 #b よつば

f1 <- function(x,y){ #x 出会ったふたば台数 #y 出会った四つ葉台数
library(gmp)
n = as.bigz(n) 
p <- factorialZ(m)*(a^x)*(b^y)*((n-a-b)^(m-x-y)) / (factorialZ(x)*factorialZ(y)*factorialZ(m-x-y)*n^m)
return(p)
}

list1 <- matrix(0, nrow = (m+1)*(m+2)/2, ncol = 2)
ct <- 1
for(j in 0:m){
      for(k in 0:m){
            if(m - j - k >= 0){
            list1[ct,] <- c(j,k)
                  ct <- ct + 1
            }
      }
}

pxy <- matrix(0, nrow = (m+1)*(m+2)/2, ncol = 3)
bigxy <- as.bigq(numeric((m+1)*(m+2)/2))
for(i in 1:((m+1)*(m+2)/2)){
      pxy[i,] <- c(list1[i,], gmp::asNumeric(f1(list1[i,1],list1[i,2])))
      bigxy[i] <- f1(list1[i,1],list1[i,2])
}


f2 <- function(x=1, y=1){
sum(bigxy[which(pxy[,3] <= pxy[which(pxy[,1]==x & pxy[,2]==y),3])])
}

return(gmp::asNumeric(f2(x,y)))
}

f0(20,1,1)
f0(30,1,1)

L2 <- 2:50
res <- numeric(length(L2))
ct1 <- 1
for(m1 in L2){
      res[ct1] <- f0(m=m1, x=1, y=1)
      ct1 <- ct1 + 1
}
plot(res)

> f0(20,1,1)
[1] 0.004473894
> f0(30,1,1)
[1] 0.009909649

20台で0.45%しか起きない珍しい事象
30台でも0.99%しか起きない珍しい事象であった。
でも30台見れば1%位の確率で起こるのね。想像より
レアじゃかなったが、まあレアだと言えよう。