読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

気になっている分布をいくつか

 

変数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="")
  }
}

f:id:High_School_Student:20130824141138j:plain
;λ= 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="")
    }
  }
}

 f:id:High_School_Student:20130824141259j:plain

y=0の発生に意図的に重みを置くことになるので、πがある程度大きい場合は、高いy=0の発生から、y=1の発生が一旦落着き、まただんだんと発生が高くなる。

y=0の発生が高い理由が既知の場合に用いる分布であろう。

 


ガンマPoisson分布(負の二項分布)

過大分散(over-dispersion)が起きているカウントデータに対する分布。
過大分散とは、分散が理論より大きい状態を指す。
Poisson分布は、理論的に平均 = 分散となるが、分散がこの仮定より大きい場合に対応するのがガンマPoisson分布と言われる。

確率密度関数は以下。

f:id:High_School_Student:20130824140940j:plain

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

f:id:High_School_Student:20130824141432j:plain 

σを小さくした場合(0.01)通常のPoisson分布に非常に近くなる。

また、σを大きくすると、Zero-infrated Poisson分布のようにy=0の発生を強調した分布になるようである。