レイドで個体値100%をとN体遭遇できる確率

ポケモンGoで伝説レイド1000回を達成したとのツイートを観測した。
期待値は単純にレイド数＊1/216で計算できるが、3体以上持っている確率とかも知りたいと思った。
Rで計算しておく。

cal.raid <- function(raid_N, p){
      res.raid <- numeric(raid_N + 1)
      
      for(i in 0:raid_N){
            res.raid[i+1] <- choose(raid_N, i) * ((1-p)^(raid_N-i)) * (p ^ i)
      }
      res.raid.cum <- rev(cumsum(rev(res.raid)))
      return(res.raid.cum)
}


for(j in 100*1:10){
      temp <- cal.raid(j, 1/216)
      barplot(temp[1:21]*100, type = "b", xlab = "100%個体ゲット数", ylab ="100%個体をその数以上持っている確率(%)",
            names.arg = 0:20 )
      title(main = paste0(c("レイドゲット数 = ", j)), collapse = "NULL")
}

例えば1000回レイドの場合、7回100%に遭遇していれば上位18.5％にはいっていることになります。
下位何％かしりたければ、知りたい遭遇数＋１を参照して100から引けばいいです。
伝説ゲット数が15体しかない自分には100なんて無縁の話ですけどね。