「DNA鑑定とバースデイ・パラドクス」にて
ns を 10^(seq(from=0,to=8,by=0.1)) に対して確率を求めているが,各 ns[i] について,毎回 1 からns[i] まで重複して計算をしている始めることになっている。無駄だ。
何を計算しているのか,よく考えよう。
> system.time({
+ # 4.7兆分の1の確率、と考えるための数字
+ q<-4.7*10^12
+
+ # 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、1万人…など)
+ ns<- round(10^(seq(from=5,to=8,by=0.1)))
+
+ birthday<-rep(0,length(ns))
+ # 対数で計算しよう
+ # 同じ型のペアが少なくとも1組は存在する確率
+ 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="すくなくとも同じ型のペアが1組ある確率")
+ })
ユーザ システム 経過
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兆分の1の確率、と考えるための数字
+ q <- 4.7e12
+
+ # 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、1万人…など)
+ ns <- round(10^(seq(from=5,to=8,by=0.1)))
+ n <- 1e8
+
+ # 対数で計算する必要はない!!
+ # 同じ型のペアが少なくとも1組は存在する確率
+ birthday2 <- (1-(cumprod((q-0:(n-1))/q)))[ns]
+ plot(log10(ns), birthday2, type="l", xlab="log10(人口)", ylab="すくなくとも同じ型のペアが1組ある確率")
+ })
ユーザ システム 経過
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
※コメント投稿者のブログIDはブログ作成者のみに通知されます