2012-07-07

すべての場合を網羅するなら計算しやすい

こちらやこちらで、バースデイ・パラドクスとか、日の確率が均等でないときへの拡張のことなどを書いている
こちらでごく少ない人数の場合をを使って式で表してみた
- $P(k=1)=1$
- $P(k=2)=1-\sum_{i=1}^n p_i^2$
- $P(k=3)=1-3\sum_{i=1}^n p_i^2 + 2 \sum_{i=1}^n p_i^3$
- $P(k=4)=1- 6\times \sum_{i=1}^n p_i^2 + 8\times \sum_{i=1}^n p_i^3 -6 \times \sum_{i=1}^n p_i^4 + 3 \times (\sum_{i=1}^n p_i^2)^2$
規則がみえなくて汚い
漸化式にしてみる
- $P=(p_i);\sum_{i=1}^n p_i=1$ をnタイプの生起確率とする
- 人のタイプがすべて異なる確率をと書くこととする
  - $T(P,k=1) = 1 ;\text{when } k=1$
  - $T(P,k=2) = 1-\sum_{i=1}^n p_i^2; \text{when } k=2$
  - - ただし $Q(-i)$ はベクトル $P$ から $p_i$ を取り除いた、長さ $n-1$ のベクトルであって、その要素の和が $1$ となるように補正したベクトル
- この式の意味は以下の通り
  - $T(P,k)$ は $T(P,k-1)$ より小さくなる。小さくなる分は、初めの $k-1$ 人はすべての異なるタイプを持つが、第 $k$ 番目の人の型が、初めの $k$ 人の型のどれかひとつに相当する場合の分である
  - この「減分」はすべての $i$ について、タイプ $i$ が２人で、残りのk-2人はタイプ $i$ 以外の $n-1$ タイプから作った $k-2$ タイプの組み合わせのすべて
  - $T(Q,k-2)$ は、タイプ $i$ を除いた $n-1$ タイプの確率の総和が1になるように補正した上で、そのような $n-1$ タイプで $k-2$ 人の型がすべて異なる組み合わせの確率をすべて足し合わせたもの(これは定義)
  - $(1-p_i)^{k-2}$ は、 $Q$ の確率を補正してあるが、「本当の確率で計算するには、 $T(Q,k-2)$ の個々の $k-2$ 要素のすべてを $(1-p_i)$ 倍する必要があることを指す
  - $p_i^2$ はある $i$ タイプについて２人が観察されている確率を考慮したもの
  - 倍するのは、こう説明できる
    - $k-1$ 人がすべて異なる型である場合の数が $n(n-1)...(n-(k-2))=\frac{n!}{(n-(k-1))!}$ 。これに対して、 $k$ 番目のタイプは、 $n$ 通り取らせることができるが、そのうち、はじめの $k-1$ 人のタイプのいずれかであるのは、 $k-1$ 通りだけ。したがって、あらたに除外されるべきは、 $(k-1)\times \frac{n!}{(n-(k-1))!}$ 通り。今、これを重複するタイプの数 $n$ ごとに取りまとめて処理する(しかもどのタイプが重複するかについては相互に対等)ので、ある重複タイプについて、 $\frac{(k-1)}{n} \times \frac{n!}{(n-(k-1))!}$ 個の要素があることになる。これらの要素は、 $n-1$ から $k-1$ 個の要素を取り出す順列のすべてに関して重複しているが、その重複の程度を求めたいのであるから、 $\frac{(k-1)}{n} \times \frac{n!}{(n-(k-1))!} / (\frac{(n-1)!}{(n-1-(k-2))!})$ が求めたい値。これを計算すると $k-1$
Rで書いてみる

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))
	}	
}

n <- 4
p <- runif(n)
p <- p/sum(p)

library(gtools)

cmb <- combinations(n,3)

P <- rep(0,length(cmb[,1]))
for(i in 1:length(cmb[,1])){
	P[i] <- prod(p[cmb[i,]])
}
sum(P)*factorial(3)
1-3*sum(p^2)+2*sum(p^3)

1-3*sum(p^2)+2*sum(p^3)-3*sum(p^2)+3*2*sum(p^3)-3*2*sum(p^4) +3*(sum(p^2))^2
1-6*sum(p^2)+8*sum(p^3)-6*sum(p^4) +3*(sum(p^2))^2

cmb <- combinations(n,4)
P <- rep(0,length(cmb[,1]))
for(i in 1:length(cmb[,1])){
	P[i] <- prod(p[cmb[i,]])
}

sum(P)*factorial(4)

my.calcP(p,1)
my.calcP(p,2)
my.calcP(p,3)
my.calcP(p,4)

当たり前ですが、すべてが等確率のときは
- $T(P=(1/n,1/n,...),k) = (\frac{1}{n})^k \frac{n!}{(n-k)!}$

q <- rep(1/n,n)
out.q <- rep(0,n)
for(i in 1:n){
	out.q[i] <-my.calcP(q,i)
}
nn <- 1:n
(1/n)^(nn) * factorial(n)/factorial(n-nn)
out.q

確率のばらつきの影響について定性的な情報を得る(ために限定したパターンによって定量する)

# 単一事象だけ高確率で残りが等確率の場合

n <- 10
p <- rep(1,n)
p <- p/sum(p)

fs <- 1:n
outs <- matrix(0,length(fs),n)
for(i in 1:length(fs)){
	f <- fs[i]
	q <- c(f,rep(1,n-1))
	q <- q/sum(q)
	for(j in 1:n){
		outs[i,j] <- my.calcP(q,j)
	}
}

matplot(t(outs),type="l")

# ２つの確率だけがあって、その比率が異なる場合

f <- 10
gs <-1:n
outs.g <- matrix(0,length(gs),n)
for(i in 1:length(gs)){
	q <- c(rep(f,gs[i]),rep(1,n-gs[i]))
	q <- q/sum(q)
	for(j in 1:n){
		outs.g[i,j] <- my.calcP(q,j)
	}
}
matplot(t(outs.g),type="l")


# 最大確率は同じで、残りの確率は等確率、ただし、事象の種類数が異なる場合

ns <- n:(n+5)
outs.n <- matrix(0,length(ns),n)

ff <- 5
p.max <- 1/n * ff
for(i in 1:length(ns)){
	q <- c(p.max,rep((1-p.max)/(ns[i]-1),ns[i]-1))
	for(j in 1:n){
		outs.n[i,j] <- my.calcP(q,j)
	}
}

matplot(t(outs.n),type="l")

さらに、すべての人が型を持つ確率として独自の確率ベクトルを持つように一般化しよう
人数と型との種類数が等しければ、非負正方行列が作れて、全員が異なる確率というのは、パーマネント(パーマネントについて)になる
さらに複雑化するとすれば、ある人が型をとる確率と別の人のそれとが独立ではないという仮定の導入