2012-06-28

すごく小さい計算

バースデイ・パラドクスというのがある(こちら)
「遺伝情報・DNA鑑定と刑事法」という文章(こちら)で、DNA鑑定に引き写して、バースデイ・パラドクスについて触れられている
この場合の確率計算は、「たまたま一致する確率が4.7兆分の１」だったり、実は、もっと小さい値で計算したいこともあって、10京分の１とか、1ガイ( $10^20$ )分の１とかのこともあるらしい
そして、一致ペアがいるかどうかの比較母集団が日本人全体(〜１億)とかだったりするらしい
とにかく確率はとても小さかったり、人数がとても多かったりする
対数を使うにしても、Rでやるにも $n!;n=10^9$ とかになってくると、メモリの問題が出てくるようだ
そんなのを少し回避しながらけいさんしてみた

# 4.7兆分の１の確率、と考えるための数字
q<-4.7*10^12

# 調べる人数(世界人口65億人だったり、日本人口1.2億人だったり、母集団100万人、10万人、１万人…など)
ns<-10^(seq(from=0,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="すくなくとも同じ型のペアが１組ある確率")

# 予備的

my.f <- function(q,n){
	sum(log(1-(1:n)/q))
}


my.f2 <- function(q,n){
	tmp <- log(1-(1:n)/q)
	cumsum(tmp)
}
n<-1*10^6
#my.f(q,n)
#tmp.out<-my.f2(q,n)

s<-sample(1:n,1000)
plot(s,tmp.out[s])


my.f3 <- function(q,n,k=8){
	kizami <- 10^k
	if(n<=kizami){
		tmp<-list()
		tmp[[1]]<-my.f2(q,n)
		return(tmp)
	}else{
		tmp<-list()
		tmp2<-list()
		resid <- n%%(kizami)
		a <- (n-resid)/(kizami)
		now <-1
		start <- 0
		for(i in 1:a){
			tmp[[i]] <- log(1-(now:(now+kizami-1))/q)
			tmp2[[i]] <- cumsum(tmp[[i]]) + start
			now <- now + kizami
			start <- tmp2[[i]][length(tmp2[[i]])]
		}
		tmp[[a+1]]<-log(1-(now:n)/q)
		tmp2[[a+1]] <- start + tmp[[a+1]]
		return(tmp2)
	}
}

#tmp.out.2 <- my.f2(10^18,10^6)
tmp.out.3 <- my.f3(10^12,10^6)

s<-sort(sample(1:n,1000))
plot(s,tmp.out.3[[1]][s])
#plot(tmp.out.2[s],tmp.out.3[[1]][s])


# 10京
q <- 10 * 10^16
n<-2*10^8
k<-5
tmp.out.1 <- my.f3(q,n,k)
kizami<-10^k
xxx1 <- yyy1 <- c()
for(i in 1:length(tmp.out.1)){
	xxx1 <- c(xxx1, s+(i-1)*kizami)
	yyy1 <- c(yyy1,tmp.out.1[[i]][s])
}

plot(xxx1,yyy1,type="l",main="10京",xlab="1.25億",ylab="自然対数")
# 日本人人数
j <- 1.25*10^8
abline(v=j,col=2)


plot(xxx1,exp(yyy1),type="l",main="10京",xlab="1.25億")
abline(v=j,col=2)


# 570京

q <- 570 * 10^16
n<-2*10^8
k<-5

tmp.out.2 <- my.f3(q,n,k)

xxx2 <- yyy2 <- c()
for(i in 1:length(tmp.out.2)){
	xxx2 <- c(xxx2, s+(i-1)*kizami)
	yyy2 <- c(yyy2,tmp.out.2[[i]][s])
}

plot(xxx2,yyy2,type="l",main="570京",xlab="1.25億",ylab="自然対数")
abline(v=1.25*10^8,col=2)

plot(xxx2,exp(yyy2),type="l",main="570京",xlab="1.25億")
abline(v=1.25*10^8,col=2)


# 10京、570京、1ガイ
q <- 10^20
n<-2*10^8
k<-5
tmp.out.3 <- my.f3(q,n,k)

xxx3 <- yyy3 <- c()
for(i in 1:length(tmp.out.3)){
	xxx3 <- c(xxx3, s+(i-1)*kizami)
	yyy3 <- c(yyy3,tmp.out.3[[i]][s])
}

plot(xxx3,yyy3,type="l",main="1ガイ",xlab="1.25億",ylab="自然対数")
# 日本人人数
j <- 1.25*10^8
abline(v=j,col=2)

plot(xxx3,exp(yyy3),type="l",main="1ガイ",xlab="1.25億")
abline(v=j,col=2)

ylim<-c(0.0,0.2)
par(mfcol=c(1,3))
plot(xxx1,1-exp(yyy1),ylim=ylim,type="l",main="10京",xlab="1.25億")
abline(v=j,col=2)
plot(xxx2,1-exp(yyy2),ylim=ylim,type="l",main="570京",xlab="1.25億")
abline(v=1.25*10^8,col=2)
plot(xxx3,1-exp(yyy3),ylim=ylim,type="l",main="1ガイ",xlab="1.25億")
abline(v=j,col=2)