ジョンとヨーコのイマジン日記

準備

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

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

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

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

$\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")

$\beta = (\beta_1, \beta_2)' = (2, 0.5)'$ としてポアソン分布でデータを生成し, $n=100$ の負の二項分布で推定しました.

折れ線グラフで見ると $\beta$ の事後分布はシミュレーションで設定した正解の近くにばらついているようです.

尤度をポアソン分布にした場合と負の二項分布にした場合で事後分布を比較します.

等高線は尤度と事前分布の積をパラメータの関数と見たものです. 定数倍の違いを無視したらこれは正確な事後分布と一致します. この場合は負の二項分布でやっても二項分布でやっても大差ないことがわかります.

点はギブスサンプリングで得られた乱数です. 事後分布に近いものが得られています.

どの程度の $n$ なら十分大きいと言えるかはデータのスケールによるでしょう.

$n$ をさらに大きくしてやると, ギブスサンプラーは自己相関が大きくなる印象です.

メトロポリス・ヘイスティングス法などで $n$ についてもデータから決めることもできるでしょう.

というか現実的には, 単純なポアソン回帰なら最尤法かMAP法でいいような気がしますが, ポアソン分布を使ったモデルを魔改造したい場合には, こういうテクニックが使える場合もあるかもしれません.

はじめに：身長の分布

身長の分布が正規分布っぽくなるというのはよく言われる。

ためしに、

学校保健統計調査－令和2年度（確定値）の結果の概要：文部科学省

から高校生男子17歳の身長の分布を持ってきて図示してみる。

保険統計調査のエクセルの表は1cm刻みで割合（千分率）になっている。

青い点線が正規分布で黒い実線は身長の分布である。

なかなか正規分布っぽい。

この青い点線の正規分布のパラメータは度数分布表から平均と標準偏差を定義通りに求め、

平均： $\bar x = \sum_i x_i p_i$

標準偏差： $s = \sum_i (x_i - \bar x)^2 p_i$

とした。ここでは階級 $x_i$ の確率（割合）を $p_i$ と表記した。

体重の分布

同じように、高校生男子17歳の体重の分布を持ってきて図示してみる。

正規分布らしくはまったくない。

通っぽい人なら仮想的なモデル人体みたいなものを考え、その身長が $s$ 倍されるとき、体積はおおむね $s^3$ に比例、体重はおおむね体積に比例する……つまり身長が正規分布なら体重は正規分布に従う確率変数を3乗した分布になるのでは！？と思うかもしれない（2乗3乗の法則 - Wikipedia）。

正規分布に従う確率変数を3乗したときの分布は求めようと思えば求められるので微分積分学などをやった人は計算してみてもいいかもしれないが、ここでは乱数を使って楽をしよう。

階級の値 $x_i$ を $p_i$ に比例する確率で10000回復元抽出して、それを1/3乗してみる。

3乗を1/3乗して、これで正規分布っぽくなる……かと思うとならない。

青い点線が正規分布で黒い実線は乱数で再現した体重の1/3乗の分布である。

たぶん人間の体は大きい人と小さい人で相似みたいにはなってないのかもしれない。

データ分析というのはこういうもので、紙の上だけで考えた理屈はこうやってあっさり退けられる。

身長の分布再考

もう一度身長の分布を見ると最頻値の171cm、170cmあたりで正規分布より少し急峻になっている。

これは偶然だろうか。

この調査の対象になった126900人だそうなので、高校3学年で割って大体40000人くらいのサンプルはいそうである。

先に求めた平均 $\bar x$ 、標準偏差 $s$ の正規分布から40000個乱数を取ってきて、1cm刻みの度数分布表を作ってみよう。

こうすることによって正規分布からのズレが偶然だとして、どれくらいの確率で起こるのかわかる。

こういったやり方をパラメトリック・ブートストラップという。

正規乱数から作った度数分布表を5000回実際の身長の分布に重ねたのが次の図である。

赤い線が実際の身長の分布、乱数から作ったほうがグレーである。

何本もあるグレーの線が赤い線より高くなることはそうそうない。

正規乱数から作った度数分布表の最大値の分布（1-経験分布関数）を次の図に示そう。

縦の点線が実際の身長の分布の171cmのところの確率で、このように値が集中するのは正規分布ではなかなかレアといえそうである。

身長は多分人間が目視で測るので170cmくらいに心理的に誘導されるみたいなことはあり得るかもしれない。

こんな具合に、正規分布に従うデータというのは自然にはまずないと思ったほうがいいと私は思う。

統計学で「正規分布で近似」というときは、多くの場合、データをパラメータ $\theta$ を持つ確率分布 $f(x;\theta)$ でモデル化したときに、 $\theta$ の推定量が正規分布に近づく、という話で、一般のデータ $x$ そのものが正規分布になるみたいなことを保証する定理はない。

この辺の区別はややこしくて、多くの人が当たる壁と思う。

R のコード

まとめて貼るぞ。

library(readxl)
library(ggplot2)
library(dplyr)
#https://www.e-stat.go.jp/stat-search/files?page=1&toukei=00400002&tstat=000001011648

dat_h <-read_excel("./Downloads/r3_hoken_tokei_02.xlsx",skip=4)
h17 <- as.numeric(dat_h$`17歳`)
h17<- ifelse(is.na(h17),0,h17)
len <- length(h17)
height <- as.numeric(gsub("cm|～","",dat_h$区分))

df_height <- data.frame(height=height[-c(1,(len-1):len)],
                        prob=h17[-c(1,(len-1):len)]/1000)

m_h <- with(df_height, sum(prob*height))
s_h <- with(df_height, sqrt(sum(prob*(height-m_h)^2)))

ggplot(df_height, aes(x=height, y=prob))+
  geom_line()+
  geom_line(data = NULL, aes(y=diff(pnorm(c(min(height)-1,height),m_h,s_h))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)
ggsave("hist_h.png")

######

dat_w <-read_excel("./Downloads/r3_hoken_tokei_03.xlsx",skip=4)
w17 <- as.numeric(dat_w$`17歳`)
w17 <- ifelse(is.na(w17),0,w17)
len <- length(w17)
weight <- as.numeric(gsub("kg|～","",dat_w$区分))
weight <- weight[-1]
df_weight <- data.frame(weight=weight[-1],
                        prob=w17[-1]/1000)
m_w <- with(df_weight, sum(prob*weight))
s_w <- with(df_weight, sqrt(sum(prob*(weight-m_w)^2)))

ggplot(df_weight, aes(x=weight, y=prob))+
  geom_line()+
  geom_line(data = NULL, aes(y=diff(pnorm(c(min(weight)-1,weight),m_w,s_w))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)
ggsave("hist_w.png")

sim_w <- sample(df_weight$weight, size = 10000, prob=df_weight$prob, replace = TRUE)


ggplot(data.frame(stat=sim_w^(1/3)), aes(stat))+
  geom_density()+
  stat_function(fun = function(x)dnorm(x, mean=mean(sim_w^(1/3)), sd=sd(sim_w^(1/3))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)+
  labs(x="stat (1/3)",y="density")
ggsave("hist_w3.png")

ggplot(data.frame(stat=sim_w^(1/2)), aes(stat))+
  geom_density()+
  stat_function(fun = function(x)dnorm(x, mean=mean(sim_w^(1/2)), sd=sd(sim_w^(1/2))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)+
  labs(x="stat (1/2)",y="density")
ggsave("hist_w2.png")

ggplot(data.frame(stat=log(sim_w)), aes(stat))+
  geom_density()+
  stat_function(fun = function(x)dnorm(x, mean=mean(log(sim_w)), sd=sd(log(sim_w))), colour="royalblue", linetype=2, linewidth=1)+
  theme_bw(16)+
  labs(x="stat (log)",y="density")
ggsave("hist_wlog.png")

######
sample_h <- function(i, n, m_h, s_h, height){
  set.seed(i)
  z <- rnorm(n,m_h,s_h)
  r <- range(z)
  b <- seq(floor(r[1]),ceiling(r[2]))
  h0 <- hist(z, breaks = b, plot = FALSE)
  data.frame(height=b[-1],prob=h0$counts/n,id=i)
}

simout <- lapply(1:5000, function(i){sample_h(i,40000, m_h, s_h, height)})
dfout <- bind_rows(simout)

modep <- group_by(dfout,id) %>% 
  summarise(maxp=max(prob))

ggplot(dfout, aes(x=height,y=prob))+
  geom_line(aes(group=id), alpha=0.002)+
  geom_line(data = df_height, colour="firebrick")+
  theme_bw()
ggsave("hist_sample.png")

ggplot(modep, aes(x=maxp))+
  geom_step(aes(y=1-after_stat(y)), stat="ecdf",pad=FALSE)+
  geom_vline(xintercept = max(prob_h), linetype=2)+
  labs(x="stat",y="prob")+
  theme_bw(16)
ggsave("ccdf_h.png")

ジョンとヨーコのイマジン日記

キョウとアンナのラヴラブダイアイリー改め、ジョンとヨーコのイマジン日記です。

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

準備

R による実装例

身長や体重の分布は正規分布か

はじめに：身長の分布

体重の分布

身長の分布再考

R のコード