ryamadaのコンピュータ・数学メモ

2012-07-09

多項式にする

バースデイ・パラドクスでは
- $P(k) = \frac{n!}{(n-k)!}\frac{1}{n^k}$ が同じ誕生日のいない確率
これを多項式化する
- $P(k+1) = P(k) \times \frac{n-k}{n}$ なので
- $P(k) = \sum_{i=0}^k a_i^{(k)} \frac{1}{n^i}$ と置くと
- $a_i^{(k+1)} = a_i^{(k)} - k\times a_{i-1}^{(k)}$ なる関係がある
これを使って計算するソースはこちら
ここの $a_i^{(k)}$ は、実は $P(k) = \frac{n}{n}\frac{n-1}{n}...\frac{n-(k-1)}{n}$ の分子のうちの $(n-1)(n-2)...(n-(k-1))$ を展開したときの係数に他ならない

my.calcP <- function(p,k){
	if(k == 1){
		return(1)
	}else if(k == 2){
		return(1-sum(p^2))
	}else{
		n <- length(p)
		tmp <- 0
		for(i in 1:n){
			q <- p[-i]
			q <- q/sum(q)
			tmp <- tmp + my.calcP(q,k-2) * (1-p[i])^(k-2) * p[i]^2 
		}
		return(my.calcP(p,k-1) - tmp * (k-1))
	}	
}
# 多項式の係数を計算する
birthday3 <- function(k){
	ret <- matrix(0,k,k)
	ret[1,1] <- 1
	for(i in 2:k){
		ret[i,] <- ret[i-1,]
		ret[i,] <- ret[i,] - (i-1) * c(0,ret[i,1:(k-1)])
	}
	ret
}

birthday3(3)

# 多項式を使ってバースデイ・パラドクス確率を計算する
birthday4 <- function(k,n){
	q <- birthday3(k)
	m <- 1/n
	M <- m^(0:(k-1))
	q %*% M
}
birthday4(5,5)
p <- rep(1,n)
p <- p/sum(p)
my.calcP(p,1)
my.calcP(p,2)
my.calcP(p,3)
my.calcP(p,4)
my.calcP(p,5)