この先、どうなる？ - ryamadaの遺伝学・遺伝統計学メモ

今、N回の観測をして、○がn1,×がn2回だったとする
この後、m回繰り返したら、m回のうち、何回成功するんだろう？
まず観測データから、○確率pの確率密度関数をベータ関数で表して $Pr(p) = \beta(n1+1,n2+1)=\frac{p^{n1}(1-p)^{n2}}{B(n1+1,n2+1)}$ 、ただし、 $B(a,b)$ はベータ関数
今、○確率がpのときにm回のうちk回成功するのは $\begin{pmatrix}m\\k\end{pmatrix} p^k (1-p)^{m-k}$ なので
結局 $\int_0^1\begin{pmatrix}m\\k\end{pmatrix} p^k (1-p)^{m-k} \times Pr(p) dp$
これを式変形すると
$\begin{pmatrix}m\\k\end{pmatrix}\frac{B(k+n1+1,m-k+n2+1)}{B(n1+1,n2+1)}\int_0^1 \frac{p^{k+n1}(1-p)^{m-k+n2}}{B(k+n1+1,m-k+n2+1)}$ となるが、積分は１なので
$\begin{pmatrix}m\\k\end{pmatrix}\frac{B(k+n1+1,m-k+n2+1)}{B(n1+1,n2+1)}$
これって、なんかの名前が付いているはずだと思うのだが…はて？

# n1の○、n2の×の下で、後m回
my.partial.binom <- function(m,n1,n2){
# ○は0,1,...,mになりえる
	M <- 0:m
	a <- lchoose(m,M)
	b <- lbeta(n1+1,n2+1)
	B <- lbeta(M+n1+1,M[length(M):1]+n2+1)
	exp(a+B-b)
}

plot(my.partial.binom(10,2,3))
plot(my.partial.binom(13,2,3))
# 既観察n1,n2
n1 <- 2
n2 <- 4
# 既観察とこれからの観察の総観察数N
N <- (n1+n2):30
# 追加試行の打ち○の回数別の確率のデータ格納
X <- matrix(0,length(N),max(N)+1)
# 別の格納庫
Ps <- list()
# ○率p
p <- seq(from=0,to=1,length=100)
# Nの値ごとに、○率pがどれくらいの確率になるかの格納
ps <- matrix(0,length(N),length(p))
for(i in 1:length(N)){
	tmp <- my.partial.binom(N[i]-(n1+n2),n1,n2)
	X[i,1:length(tmp)] <- tmp
	Ps[[i]] <- tmp
	for(j in 1:length(tmp)){
		ps[i,] <- ps[i,] + Ps[[i]][j] * dbeta(p,n1+(j-1)+1,n2+N[i]-(n1+n2)-(j-1)+1)
	}
	#print(sum(tmp*(0:(length(tmp)-1))))/(N-(n1+n2))
}

matplot(0:max(N),t(X),type="l")
# Nが変わっても○率の分布は変わらない
matplot(p,t(ps),type="l")
# Nが変わっても、○率の期待値は変わらない
for(i in 1:length(N)){
	print(sum(ps[i,] * p))
}