ポリア・ガンマ分布と負の二項分布を用いたポアソン回帰の近似的なギブスサンプラー

準備

多項ロジスティック回帰について

ポリア・ガンマ分布を用いた多項ロジスティック回帰のギブスサンプラー - ジョンとヨーコのイマジン日記

に書いた. たぶんこちらを先に読んだほうがいい.

ここではポアソン回帰, つまり次のモデルを考える.

$\displaystyle y_i \sim \mathrm{Pois}(\exp(x_i' \beta))$

ところで, ポアソン分布は平均と分散が同じなので, ポアソン分布より分散の大きいカウントデータのモデリングには負の二項分布が使われることがある.

負の二項分布は次の確率関数を持つ.

$\displaystyle \mathrm{NegBin}(y; n, \mu)= \frac{\Gamma(y+n)}{\Gamma(n)\Gamma(y+1)}\left(\frac{\mu}{n+\mu}\right)^n \left(\frac{n}{n+\mu}\right)^y$

（パラメータの置き方は統計ソフトによって微妙に違うことがある.）

負の二項分布は $n \to \infty$ の極限でポアソン分布になる. たとえば,

mathematical statistics - Poisson as a limiting case of negative binomial - Cross Validated

を参照してほしい.

ここから, $n$ に適当に大きい数字を入れて負の二項分布の尤度を計算すれば, ポアソン分布と大体一致するのではないかと思いつく.

この方針で少し負の二項分布の式を変形してみる.

$\mu = \exp(\psi - \log n)$ となるように $\psi$ を置くと,

$\displaystyle \mu/(\mu+n) = \exp(\psi - \log n) / (1 + \exp(\psi - \log n))$

さらに $\eta = \psi - log n$ と改めて置くと, 負の二項分布のサンプル1つあたりの尤度 $L$ は,

$\displaystyle L \propto \frac{\exp(\eta)^{y}}{\{1+\exp(\eta)\}^{y+n}}$

と表せ, 結果, ポリア・ガンマ分布の性質がそのまま使える.

$\psi = x_i' \beta$ の形で回帰構造を入れ, $\beta_j$ の事前分布を正規分布

$\displaystyle \beta \sim \mathcal{N}(0, \tau^{-1})$

とすると, このとき full conditional はポリア・ガンマ分布と多変量正規分布を用いて、

$\displaystyle \omega_{i} \mid \beta \sim \mathrm{PG}(y_i+n, x_i \beta_j)$

$\displaystyle \beta \mid \omega_i \sim \mathcal{N}(\hat \mu, \hat V)$

となる. ここで,

$\displaystyle \hat V = (X'\Omega X+ \tau )^{-1}$

$\displaystyle \hat \mu = \hat V(X'(y-0.5(y+n)+\Omega \log n))$

と置いた.

また, $y = (y_1, \ldots ,y_N)'$ , $\displaystyle \Omega$ は $\omega_{i}$ を対角成分とする対角行列とした.

R による実装例

library(BayesLogit)
library(ggplot2)

gibbs_nbinom <- function(y,n,X,tau,iter){
  N <- length(y)
  d <- ncol(X)
  B <- diag(tau, d)
  beta_tilde <- matrix(0, iter, d)
  beta_tilde[1,] <- rnorm(d)
  logn <- log(n)
  m <- y + n
  for(i in 2:iter){
    kappa <- 0.5*y-0.5*exp(logn)
    eta <- X%*%beta_tilde[i-1,]-logn
    omega <- rpg(N, m, eta)
    Omega <- diag(omega)
    Vinv <- t(X)%*%Omega%*%X + B
    U <- chol(Vinv)
    A <- forwardsolve(t(U), t(X)%*%(kappa+logn*omega))
    mu <- backsolve(U, A) #equivalent to #mu <- solve(Vinv%*%(t(X)%*%(y-0.5*n)))
    beta_tilde[i,]  = mu + backsolve(U, rnorm(d)) #multiply to inverse of U
  }
  return(list(beta=beta_tilde))
}

set.seed(1234)
beta <- c(2,0.5)
X <- matrix(runif(2*100),100,2)
y <- rpois(100,exp(X%*%beta))
print(max(y))
out_sample <- gibbs_nbinom(y=y,n=100,X=X,tau=1,iter=4000)

dfs <- expand.grid(iter=1:4000, index=1:2)
dfs$value <- as.vector(out_sample$beta)

dft <- data.frame(beta=beta, index=1:2)
ggplot(dfs, aes(x=iter, y=value))+
  geom_line(colour="darkgrey")+
  geom_hline(data=dft, aes(yintercept=beta), colour="royalblue")+
  facet_grid(index~., scales = "free", labeller = label_both)+
  theme_classic(14)+
  theme(strip.text.y = element_text(angle = 0),
        axis.text = element_text(colour = "black"))

#ggsave("tracep.png")

lp_pois <- function(beta){
  sum(dpois(y,exp(X%*%beta), log = TRUE))+sum(dnorm(beta,log=TRUE))  
}

lp_negbin <- function(beta,n){
  sum(dnbinom(y,size=n,mu=exp(X%*%beta), log = TRUE))+sum(dnorm(beta,log=TRUE))  
}


burnin <- 1:500

dfp <- data.frame(out_sample$beta[-burnin,])
colnames(dfp) <- c("b1","b2")

df <- expand.grid(b1=seq(1.5,2.5,by=0.005),b2=seq(0,1,by=0.005))
lp_v_p <- apply(df, 1, function(b){lp_pois(c(b[1],b[2]))})
lp_v_nb <- apply(df, 1, function(b){lp_negbin(c(b[1],b[2]), 100)})

ggplot(df, aes(x=b1, y=b2))+
  geom_contour(aes(z=exp(lp_v_p), colour="poisson"))+
  geom_contour(aes(z=exp(lp_v_nb), colour="negbin"))+
  geom_point(data=dfp, alpha=0.1)+
  scale_colour_manual(values = c("royalblue","orange2"))+
  labs(colour="likelihood")+
  theme_classic(16)

#ggsave("post.png")