廿TT

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

オフセット項のある負の二項分布

今日の川柳

モデル

パラメータがガンマ分布に従って変化するようなポアソン分布は負の二項分布になることが知られています。

オフセット(既知であり推定しない変数)\tau のあるポアソン分布のパラメータがガンマ分布に従って変化する場合、どのような分布になるかを考えます。

 y \sim \mathrm{Poisson}(\lambda \tau)
 \lambda \sim \mathrm{Gamma}(\alpha, \beta)

\lambda積分消去します。

 \displaystyle \int^{\infty}_{0} \frac{(\lambda \tau)^y}{y!}e^{-\lambda \tau} \times \frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda^{\alpha-1}e^{-\beta\lambda}\,d\lambda

\lambda に依存しない因子を積分記号の外に出し、r=(\tau+\beta)\lambda とおくと、

 \displaystyle p(y)=\frac{\tau^y\beta^\alpha}{y!\Gamma(\alpha)}\left(\frac{1}{\tau+\beta}\right)^{y+\alpha}\int^{\infty}_{0} r^{y+\alpha}e^{-r}\,dr\\
=\displaystyle \frac{\Gamma(y+\alpha)}{y!\Gamma(\alpha)}\left(\frac{1}{\tau+\beta}\right)^{y+\alpha}\tau^y\beta^\alpha

 \tau = 1 のとき、

 \displaystyle p(y)=\displaystyle \frac{\Gamma(y+\alpha)}{y!\Gamma(\alpha)}\left(\frac{\beta}{\beta+1}\right)^\alpha\left(\frac{1}{\beta+1}\right)^y

となり、普通の負の二項分布と一致します。

ケーススタディ

統計局ホームページ/第六十三回日本統計年鑑 平成26年−第2章 人口・世帯 の「2 - 5 都道府県別人口集中地区人口,面積及び人口密度(エクセル:34KB)」を使い、面積をオフセットとして人工の分布を調べます。

まずはオフセットのあるポアソン分布を用いた場合。

最尤推定されたパラメータを持つ分布から1万回乱数の抽出をおこない、95%区間をプロットしてみます。

縦軸が人工、横軸が面積です。

f:id:abrahamcow:20180806053935p:plain

平均はだいたい捉えていますが、データの散らばり具合をカバーしきれていないことがわかります。

次にオフセットのある負の二項分布を用いた場合。

f:id:abrahamcow:20180806054135p:plain

データの散らばりに近い分布が推定されていることがわかります。

R のコード

library(tidyverse)
library(readxl)
dat <-read_xls("~/Downloads/y0205000.xls",skip=8) %>% 
  slice(-c(1:8)) %>% 
  select(X__1,人口,面積) %>% 
  set_names(c("pref","pop","area")) %>% 
  dplyr::filter(!is.na(pop)) %>% 
  mutate(pop=as.integer(pop),area=as.numeric(area))

fit1 <- glm(pop~1,offset=log(area),family=poisson,data = dat)
sim1 <-matrix(rpois(10000*47,dat$area*exp(fit1$coefficients)),47,10000)
int1 <-as.data.frame(t(apply(sim1,1,quantile,prob=c(0.025,0.975)))) %>% 
  set_names(c("lower","upper")) %>% 
  mutate(area=dat$area)

ggplot(dat,aes(x=area))+
  geom_point(aes(y=pop),pch=1)+
  geom_ribbon(data=int1,aes(ymin=lower,ymax=upper),alpha=0.3)+
  theme_bw()

dnbinom2 <- function(y,alpha,beta,tau=1,log=FALSE){
  lp <-lgamma(y+alpha) - lfactorial(y) - lgamma(alpha) - 
    (y+alpha)*log(tau+beta)+
    y*log(tau) + alpha*log(beta)
  if(log){
    return(lp)
  }else{
    return(exp(lp))
  }
}
ll <- function(par,y,tau){
  sum(dnbinom2(y,exp(par[1]),exp(par[2]),tau,log = TRUE))
}

opt <-optim(c(0,0),ll,y=dat$pop,tau=dat$area,control = list(fnscale=-1))

lambda <- rgamma(10000*47,exp(opt$par)[1],exp(opt$par)[2])
sim2 <-matrix(rpois(10000*47,dat$area*lambda),47,10000)
int2 <-as.data.frame(t(apply(sim2,1,quantile,prob=c(0.025,0.975)))) %>% 
  set_names(c("lower","upper")) %>% 
  mutate(area=dat$area)

ggplot(dat,aes(x=area))+
  geom_point(aes(y=pop),pch=1)+
  geom_ribbon(data=int2,aes(ymin=lower,ymax=upper),alpha=0.3)+
  theme_bw()

参考文献

オフセットや負の二項分布については『データ解析のための統計モデリング入門』に解説があります。

追加情報