ダメ出し：ベクトル計算をしよう

「DNA鑑定とバースデイ・パラドクス」にて

ns を 10^(seq(from=0,to=8,by=0.1)) に対して確率を求めているが，各 ns[i] について，毎回 1 からns[i] まで重複して計算をしている始めることになっている。無駄だ。
何を計算しているのか，よく考えよう。

> system.time({
+ # 4.7兆分の１の確率、と考えるための数字
+ q<-4.7*10^12
+
+ # 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、１万人…など)
+ ns<- round(10^(seq(from=5,to=8,by=0.1)))
+
+ birthday<-rep(0,length(ns))
+ # 対数で計算しよう
+ # 同じ型のペアが少なくとも１組は存在する確率
+ for(i in 1:length(ns)){
+     birthday[i]<-1-exp(sum(log(1-1:(ns[i]-1)/q)))
+ }
+
+ plot(log10(ns),birthday,type="l",xlab="log10(人口)",ylab="すくなくとも同じ型のペアが１組ある確率")
+ })
   ユーザ   システム       経過  
    22.124     10.256     32.140

i = 1～1e8 まで，「1刻み」で各項に付加される確率 (q-i)/q, i=0, 1, 2, ... を計算しておき，cumprod して余事象の確率を求め，1 から引く。
さらに，この場合は対数で計算する必要はさらさらない（大きな数は出てこない）。
以下のようにすれば，計算時間は 1/3 になる。

> system.time({
+ # 4.7兆分の１の確率、と考えるための数字
+ q <- 4.7e12
+
+ # 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、１万人…など)
+ ns <- round(10^(seq(from=5,to=8,by=0.1)))
+ n <- 1e8
+
+ # 対数で計算する必要はない！！
+ # 同じ型のペアが少なくとも１組は存在する確率
+ birthday2 <- (1-(cumprod((q-0:(n-1))/q)))[ns]
+ plot(log10(ns), birthday2, type="l", xlab="log10(人口)", ylab="すくなくとも同じ型のペアが１組ある確率")
+ })
   ユーザ   システム       経過  
     6.135      3.136      9.195

> all.equal(birthday, birthday2)
[1] TRUE

> cbind(birthday, birthday2, birthday-birthday2)
         birthday   birthday2              
 [1,] 0.001063254 0.001063254  0.000000e+00
 [2,] 0.001684635 0.001684635  0.000000e+00
 [3,] 0.002668625 0.002668625  0.000000e+00
 [4,] 0.004226196 0.004226196  0.000000e+00
 [5,] 0.006689827 0.006689827  0.000000e+00
 [6,] 0.010581894 0.010581894  0.000000e+00
 [7,] 0.016719167 0.016719167  0.000000e+00
 [8,] 0.026368242 0.026368242 -1.110223e-16
 [9,] 0.041467409 0.041467409  0.000000e+00
[10,] 0.064919822 0.064919822  1.110223e-16
[11,] 0.100919658 0.100919658  0.000000e+00
[12,] 0.155157813 0.155157813  0.000000e+00
[13,] 0.234496703 0.234496703  0.000000e+00
[14,] 0.345260597 0.345260597  0.000000e+00
[15,] 0.488920866 0.488920866  0.000000e+00
[16,] 0.654868549 0.654868549  0.000000e+00
[17,] 0.814751459 0.814751459  0.000000e+00
[18,] 0.930901321 0.930901321  0.000000e+00
[19,] 0.985522844 0.985522844  0.000000e+00
[20,] 0.998784153 0.998784153  0.000000e+00
[21,] 0.999976020 0.999976020  0.000000e+00
[22,] 0.999999952 0.999999952  0.000000e+00
[23,] 1.000000000 1.000000000  0.000000e+00
[24,] 1.000000000 1.000000000  0.000000e+00
[25,] 1.000000000 1.000000000  0.000000e+00
[26,] 1.000000000 1.000000000  0.000000e+00
[27,] 1.000000000 1.000000000  0.000000e+00
[28,] 1.000000000 1.000000000  0.000000e+00
[29,] 1.000000000 1.000000000  0.000000e+00
[30,] 1.000000000 1.000000000  0.000000e+00
[31,] 1.000000000 1.000000000  0.000000e+00