メモリと計算時間のトレードオフ

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

・対数を使うにしても、Rでやるにもn!;n=10^9とかになってくると、メモリの問題が出てくるようだ
・そんなのを少し回避しながら計算してみた

メモリと計算時間は互いにトレードオフの関係にある。
計算時間は掛かる（べらぼうに！でも許容範囲）がメモリは全く使わないプログラムを書けばよいだけだ。
このプログラムは引用元のプログラムが以下に無駄な計算をしているかを明示している。

func <- function(n) {
    q <- 4.7e12
    p <- 1
    for (i in 1:n - 1) {
        p <- p*(q-i)/q
    }
}

n = 10^i で i が 1 大きくなると，計算時間は 10 倍になる。
システム時間がほとんど掛かっていないところにも注目。

> for (i in 0:8) {
+     print(system.time(func(10^i)))
+ }
   ユーザ   システム       経過  
     0.001      0.000      0.000
   ユーザ   システム       経過  
         0          0          0
   ユーザ   システム       経過  
     0.001      0.000      0.001
   ユーザ   システム       経過  
     0.002      0.000      0.002
   ユーザ   システム       経過  
     0.018      0.000      0.018
   ユーザ   システム       経過  
     0.190      0.001      0.190
   ユーザ   システム       経過  
     1.941      0.011      1.939
   ユーザ   システム       経過  
    19.492      0.089     19.387
   ユーザ   システム       経過  
   196.289      1.623    196.516

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

「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