気になっている分布をいくつか
変数yがカウントデータ（y = 0, 1, 2, 3, ….）のときの分布に関していくつか。
Poisson分布
http://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%88%86%E5%B8%83
パラメータλのみによって形状が定まる。
平均と分散は等しい
E(y) =λ
Var(y) =λ
λ=1,2,3,4,5と変化させつつ、y = 0,1,2,3,…,10のときの発生率をプロット。
lam <- c(1,2,3,4,5)
y <- 0:10
# 行列に値を格納する。行がλ=1,2,3,4,5、列がy=0,1,2,...,10。
Res <- 1:(length(lam)*length(y))
dim(Res) <- c(length(lam), length(y))
for(i in 1:length(lam)){
for(j in 1:length(y)){
# dpois()はPoisson分布の確率密度を返す関数。
Res[i,j] <- dpois(y[j], lambda=lam[i])
}
}
# 重ね合わせプロット。λ=1のときのプロットにλ=2,3,4,5を重ねるのでif()を使用。
for(i in 1:length(lam)){
if(i == 1){
plot(y, Res[i,], type="o", col=i, ylim=c(0,0.45), main="Poisson分布, λ=1,2,3,4,5",xlab="y", ylab="prob(y)" )
}else{
par(new=TRUE)
plot(y, Res[i,], type="o", col=i, ylim=c(0,0.45), xlab="", ylab="")
}
}
黒；λ= 1 赤：λ= 2 緑：λ= 3 青：λ= 4 水色：λ= 5
Zero-infrated Poisson分布
http://en.wikipedia.org/wiki/Zero-inflated_model
0カウントの発生率を多めに見積もる分布。
Prob(y) = π+ (1-π) f( y ) if y = 0
Prob(y) = (1-π) f( y ) if y = 1, 2, 3, 4,…
f( y )は、Poisson分布などの確率密度分布を示す。
πは、上リンクでは“probability of extra zeros”とされており、y=0の発生率がy>0の場合に比べπだけ高くする、という役割を果たすパラメータ。0から1の間をとり、π= 0とすると、y = 0, 1, 2, 3, 4,… において単純にProb(y) = f( y )となる。
Zero-infrated Poisson分布も上リンクに記載されている。
E(y) =λ(1 - π)
Var(y) =λ(1 - π)(1 + πλ)
と定義される。
Zero-infrated Poisson分布をπ=0.1, 0.2, 0.4. 0.7別にプロットしてみる。
pi1 <- c(0.1, 0.2, 0.4, 0.7)
lam <- c(1,2,3,4,5)
y <- 0:10
# 行列に値を格納する。行がλ=1,2,3,4,5、列がy=0,1,2,...,10。
Res <- 1:(length(lam)*length(y))
dim(Res) <- c(length(lam), length(y))
# さらに、π=0.1,0.2,0.4,0.7別に行列をリストで管理。
for(i in 1:length(pi1)){
if(i == 1){
LRes <- list(Res)
}else{
# append()でリストに要素を追加
LRes <- append(LRes, list(Res))
}
}
# ZI Poisson分布の確率密度を取得。π=0.1,0.2,0.4,0.7のときのそれぞれの値がリストの要素。
for(k in 1:length(pi1)){
for(i in 1:length(lam)){
for(j in 1:length(y)){
# y=0のとき、πだけ発生率が水増しされる。
if(y[j] == 0){
LRes[ [k] ][i,j] <- pi1[k] + (1-pi1[k]) * dpois(y[j], lambda=lam[i])
}else{
LRes[ [k] ][i,j] <- (1-pi1[k]) * dpois(y[j], lambda=lam[i])
}
}
}
}
# π別にプロットを描く。
par(mfrow=c(2,2))
for(k in 1:length(pi1)){
# グラフタイトルを別で作っておく。
title <- paste("ZI Poisson分布 π=", pi1[k] , ", λ=1,2,3,4,5")
for(i in 1:length(lam)){
if(i == 1){
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.9), main=title,xlab="y", ylab="prob(y)" )
}else{
par(new=TRUE)
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.9), xlab="", ylab="")
}
}
}
y=0の発生に意図的に重みを置くことになるので、πがある程度大きい場合は、高いy=0の発生から、y=1の発生が一旦落着き、まただんだんと発生が高くなる。
y=0の発生が高い理由が既知の場合に用いる分布であろう。
ガンマPoisson分布（負の二項分布）
過大分散（over-dispersion）が起きているカウントデータに対する分布。
過大分散とは、分散が理論より大きい状態を指す。
Poisson分布は、理論的に平均 = 分散となるが、分散がこの仮定より大きい場合に対応するのがガンマPoisson分布と言われる。
確率密度関数は以下。
E(y) = μ
Var(y) = μ + σ (μ^2)
平均と分散が、σ（過大分散パラメータ）の指定によって独立している。
http://ja.wikipedia.org/wiki/%E8%B2%A0%E3%81%AE%E4%BA%8C%E9%A0%85%E5%88%86%E5%B8%83
負の二項分布の確率関数を、パラメータを置き換え、Combinationを階乗で表現し、階乗をガンマ関数に置き換えると上のガンマPoissonに変形できる。
ガンマPoisson分布をσ=0.01, 1, 2, 4別にプロットしてみる。
mu <- c(1, 2, 3, 4, 5)
sig <- c(0.01, 1, 2, 4)
y <- 0:10
# 行列に値を格納する。行がμ=1,2,3,4,5、列がy=0,1,2,...,10。
Res <- 1:(length(mu)*length(y))
dim(Res) <- c(length(mu), length(y))
# さらに、σ=0.01,1,2,4別に行列をリストで管理。
for(i in 1:length(sig)){
if(i == 1){
LRes <- list(Res)
}else{
LRes <- append(LRes, list(Res))
}
}
# Gamma Poissonを関数（GP）を定義しておく。μ、σ、yが引数で、確率密度を返す。標準関数で準備されているかは不明。
GP <- function(mu1, sig1, y1){
f <- ( gamma(y1+(1/sig1)) / (gamma(y1+1)*gamma(1/sig1)) )*( (mu1*sig1)^y1 / (1+mu1*sig1)^(y1+(1/sig1)) )
return(f)
}
# ガンマPoisson分布の確率密度を取得。σ=0.01,1,2,4のときのそれぞれの値がリストの要素。
for(k in 1:length(sig)){
for(i in 1:length(mu)){
for(j in 1:length(y)){
LRes[ [k] ][i,j] <- GP(mu[i], sig[k], y[j])
}
}
}
# σ別にプロットを描く。
par(mfrow=c(2,2))
for(k in 1:length(sig)){
title <- paste("Γ Poisson分布 σ=", sig[k] , ", μ=1,2,3,4,5")
for(i in 1:length(mu)){
if(i == 1){
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.8), main=title,xlab="y", ylab="prob(y)" )
}else{
par(new=TRUE)
plot(y, LRes[ [k] ][i,], type="o", col=i, ylim=c(0,0.8), xlab="", ylab="")
}
}
}
σを小さくした場合（0.01）通常のPoisson分布に非常に近くなる。
また、σを大きくすると、Zero-infrated Poisson分布のようにy=0の発生を強調した分布になるようである。
