廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

ガンマ・ポアソン分布で推定しても負の二項分布で推定しても結果は一緒になるか

今日の川柳

問:ガンマ・ポアソン分布で推定しても負の二項分布で推定しても結果は一緒になるか
答:なりそう


以下のデータ生成過程を考える。
 z \sim \mathrm{Gamma}(\alpha,\alpha)
 y \sim \mathrm{Poisson}(z\lambda)

z を積分消去すると、y は負の二項分布に従うことがわかる(証明略)。

また、z の事後分布は以下のガンマ分布になる。
\mathrm{Gamma}(y+\alpha,\alpha+\lambda)(証明略)。

パラメータ z と α、λ を推定する Stan のコードはこう。

//nbinom.stan
data{
  int N;
  int y[N];
}
parameters{
  real<lower=0> alpha;
  real<lower=0> lambda;
  vector<lower=0>[N] z;
}
model{
  z ~ gamma(alpha,alpha);
  y ~ poisson(z*lambda);
}

z を積分消去してパラメータ α と λ を推定する Stan のコードはこう。

//nbinom2.stan
data{
  int N;
  int y[N];
}
parameters{
  real<lower=0> alpha;
  real<lower=0> lambda;
}
model{
  y ~ neg_binomial_2(lambda,alpha);
}
generated quantities{
  real z[N];
  for(i in 1:N){
    z[i] = gamma_rng(y[i]+alpha,alpha+lambda);
  }
}

nbinom.stan による z の分布を赤、nbinom2.stan でサンプリングした z の分布を青で書いて重ねた。

nbinom2.stan のほうが若干はやいです。

f:id:abrahamcow:20180804000326p:plain

両者はかなり一致し、赤と青が混ざって赤紫っぽくなってる。

パラメータ α と λ の事後分布も近いものが求まりました。

f:id:abrahamcow:20180804000848p:plain

f:id:abrahamcow:20180804000834p:plain

R のコードを貼ります。

require(rstan)
require(tidyverse)
rstan_options(auto_write = TRUE)
set.seed(1234)
N <- 49
z <- rgamma(N,2,2)
y <- rpois(N,z*10)
dat4stan <- list(N=N,y=y)
nbin1 <-stan_model(file="nbinom.stan")
nbin2 <-stan_model(file="nbinom2.stan")
fit1 <-sampling(nbin1,dat4stan)
fit2 <-sampling(nbin2,dat4stan)
#traceplot(fit1)
#traceplot(fit2)
ex1 <- rstan::extract(fit1)
ex2 <- rstan::extract(fit2)
z1df <-as.data.frame(ex1$z) %>% 
  gather()
z2df <-as.data.frame(ex2$z) %>% 
  gather()
ggplot(z1df,aes(x=value))+
  geom_histogram(bins=40,fill="red",alpha=0.5)+
  geom_histogram(data=z2df,bins=40,fill="blue",alpha=0.5)+
  facet_wrap(~key,scales = "free")+
  theme(strip.background = element_blank(),strip.text = element_blank(),
        panel.background = element_blank())

plot(density(ex2$alpha),col="blue",main="alpha")
lines(density(ex1$alpha),col="red")

plot(density(ex1$lambda),col="red",main="lambda")
lines(density(ex2$lambda),col="blue")