データを見て予測する - ryamadaの遺伝学・遺伝統計学メモ

こちらでポアソン分布を気にしている
その続き
L回の観察機会のうち、 $k$ ([tex:k<
この事象をN回観察する機会があったら、何回起きると予想するのだろうか？
L回の観察、k回の生起、から、ポアソン分布のパラメタ $\lambda$ の確率密度分布を、ガンマ分布として推定することにする
ガンマ分布を選ぶのは、ポアソン分布の共役事前分布だから

L <- 50
k <- 3
lambda <- seq(from=0,to=1,length=1000)
d.lambda <- dgamma(lambda,k,L)
plot(lambda,d.lambda,type="l")

いろいろな $\lambda$ について、N回の観察をしたときにt回生起する確率はポアソン分布に従う

N <- 200
t <- 0:20
seiki <- matrix(0,length(lambda),length(t))
for(i in 1:length(lambda)){
	seiki[i,] <- dpois(t,lambda[i])
}
apply(seiki,1,sum) # すべて1になる
matplot(t(seiki),type="l")
matplot(t(seiki[sample(1:length(lambda),10),]),type="l")

$\lambda$ の値を適当に細かくとって、それぞれの $\lambda$ の値に上で求めた、確率の重みをつければ、N回の観察でt回の生起が起きる確率を求めることができる

x <- matrix(0,length(lambda),length(t))
for(i in 1:length(lambda)){
	x[i,] <- d.lambda[i] * seiki[i,] * (lambda[2]-lambda[1])
}

sum(x)
seiki2 <- apply(x,2,sum)

plot(t,seiki2,type="b")

$0 \le \lambda \le \infty$ の値を連続的にとって、その全部について $\lambda$ の確率に応じて積分しよう
$\lambda$ の確率密度分布を[tex:G(k+a,L+b) = \frac{\lambda^{(k+a)-1}exp(-(L+b)\lambda)}{\Gamma*1(\frac{1}{(L+b)})^k}]と推定することにする
$\lambda$ のポアソン分布でN回の観察をすれば、 $N\lambda$ 回の生起が期待されるから、t回生起する確率は $P(N,\lambda) = \frac{(N\lambda)^t}{t!}exp(-N\lambda)$
両方の分布を $\lambda$ について積分すればよいから
$\int_{0}^{\infty} G(k+a,L+b)P(N,\lambda) d\lambda$ が知りたい確率
[tex:G(k+a,L+b)P(N,\lambda) = \frac{\lambda^{(k+a)-1}exp(-(L+b)\lambda)}{\Gamma*2(\frac{1}{(L+b)})^k}\frac{(N\lambda)^t}{t!}exp(-N\lambda)]
[tex:\frac{N^t}{\Gamma*3(\frac{1}{L+b})^k t!} \lambda^{k+a+t-1} exp(-(L+b+N)\lambda)]
ここでについて考える
- ガンマ分布の確率密度関数が $\frac{x^{t-1}exp(-\frac{x}{s})}{\Gamma(t)s^k$ であることから、 $\int_{0}^{\infty} \frac{x^{t-1}exp(-\frac{x}{s})}{\Gamma(t)s^k} = 1$
- したがって $\int_{0}^{\infty} x^{t-1}exp(-sx) dx = \Gamma(t)s^t$
これを用いると
$\int_0^{\infty} G(k+a,L+b)P(N,\lambda) d\lambda$
[tex:=\frac{N^t}{\Gamma*4(\frac{1}{L+b})^k t!} \Gamma(k+a+t) (\frac{1}{L+b})^{k+a+t}]

my.fromgammatopois <- function(N,ts,k,theta){
	ps <- (k+ts-1)
	q <- -(N+theta)
	exp(ts*log(N)-lgamma(k)-k*log(1/theta)-lgamma(ts+1) + lgamma(ps+1) - (ps+1)*log(-q))

}
a <- 1
b <- 1
k <- 2
L <- 39
N <- 100
t <- 0:100
ooo <- my.fromgammatopois(N,t,k+a,L+b)
plot(ooo,type="h")
sum(ooo)

*1:k+a

*2:k+a

*3:k+a

*4:k+a