- 今、N回の観測をして、○がn1,×がn2回だったとする
- この後、m回繰り返したら、m回のうち、何回成功するんだろう?
- まず観測データから、○確率pの確率密度関数をベータ関数で表して
、ただし、
はベータ関数
- 今、○確率がpのときにm回のうちk回成功するのは
なので
- 結局
![\int_0^1\begin{pmatrix}m\\k\end{pmatrix} p^k (1-p)^{m-k} \times Pr(p) dp](https://chart.apis.google.com/chart?cht=tx&chl=%5Cint_0%5E1%5Cbegin%7Bpmatrix%7Dm%5C%5Ck%5Cend%7Bpmatrix%7D%20p%5Ek%20%281-p%29%5E%7Bm-k%7D%20%5Ctimes%20Pr%28p%29%20dp)
- これを式変形すると
となるが、積分は1なので
![\begin{pmatrix}m\\k\end{pmatrix}\frac{B(k+n1+1,m-k+n2+1)}{B(n1+1,n2+1)}](https://chart.apis.google.com/chart?cht=tx&chl=%5Cbegin%7Bpmatrix%7Dm%5C%5Ck%5Cend%7Bpmatrix%7D%5Cfrac%7BB%28k%2Bn1%2B1%2Cm-k%2Bn2%2B1%29%7D%7BB%28n1%2B1%2Cn2%2B1%29%7D)
- これって、なんかの名前が付いているはずだと思うのだが…はて?
![f:id:ryamada22:20140623153624j:image f:id:ryamada22:20140623153624j:image](https://cdn-ak.f.st-hatena.com/images/fotolife/r/ryamada22/20140623/20140623153624.jpg)
my.partial.binom <- function(m,n1,n2){
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 <- 2
n2 <- 4
N <- (n1+n2):30
X <- matrix(0,length(N),max(N)+1)
Ps <- list()
p <- seq(from=0,to=1,length=100)
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)
}
}
matplot(0:max(N),t(X),type="l")
matplot(p,t(ps),type="l")
for(i in 1:length(N)){
print(sum(ps[i,] * p))
}