ポアソン分布の2つの起源；島谷（2017）から

趣旨：島谷健一郎『ポアソン分布・ポアソン回帰・ポアソン過程』（近代科学社；以下ではポアソン³と略記する）は、最初「トピックを絞ったおもしろい読みもの」的な内容だと思っていたが、読み直すたびに評価が上がってきて、全体的に統計学入門にいいのではないかと思うまでに至ったので、勝手に応援します。

ポアソン分布・ポアソン回帰・ポアソン過程 統計スポットライト・シリーズ

さて、ポアソン³の第1章ではこんな感じのシミュレーションが紹介される。

このGIFアニメがなにをやっているかというと、大きい空間の中にランダムに（一様分布の直積で）点を配置して、中央の四角い区間の中に入った点の数を数えている。右のパネルが点の配置で、左のパネルが集計した結果の棒グラフ。左のグレーの点線はポアソン分布の確率関数で、棒グラフがこれに近い形になる。

もう一つ：

このGIFアニメは、1か0かがランダムに（ベルヌーイ分布で）決まる格子を並べて、全体の1の数を数えている。右のパネルが格子で左のパネルが集計した結果の棒グラフ。左のグレーの点線はポアソン分布の確率関数で、棒グラフがこれに近い形になる。

これらは要は「二項分布のポアソン近似」をやっている。

こんな感じでシミュレーションや具体例を通じてポアソン分布・ポアソン回帰・ポアソン過程を理解しようとするところがポアソン³の特色である。

島谷は「本書のようなシミュレーションは『使うためだけ』の理解と『数学としての理解』の中間に位置する」と述べる（p.95; 4章の最後の方。ちなみに4章もカルバック・ライブラー情報量の直感的な意味という非常に大切なことが述べられている章だ）。私はこれは控えめに過ぎる表現だと思う。

「二項分布のポアソン近似」を（統計学の授業などで）紙とペンだけで勉強して実際にイメージを持つことは難しい。

「ポアソンの小数の法則だから数が小さいときじゃないとポアソン分布使っちゃいけないのかなー」と思ってしまう場合すらあるだろう。

シミュレーションをやってみた経験があれば「でも格子の数をすごく増やせば…」というイメージを持ちやすくなる。

また、統計学の定理として数学的に証明されるような結果はサンプルサイズ の場合の結果であったりする。

データがいっぱいあればその分だけ対象が正確にわかるようになってほしいので、もちろんそれらの結果は大事なものなのだが（こういうのを漸近論といったりする）、現実の有限の標本でどのような推定が得られるかとかは手計算するのがとても難しく、シミュレーションに頼らざるを得ない場合が多い。

つまりこういったシミュレーションを作れるテクニックを習得しておくと、ポアソン分布・ポアソン回帰・ポアソン過程に限らず、統計学全体で便利だと思う。

もっというと、なんかすごい発想とかがなくても「愚直にいろんなシミュレーションをやってみました」が、普通に現実で役に立つ研究であることもある。

また、「この分布からこの変数が決まって、次に…」という「データ生成過程」に対する素朴な理解は統計モデリングの際に重要だと思うが、モデルでデータをまねっこするということを明示的に教えてもらえる場面は、実は少ないのじゃないかとも思う。

ポアソン³は主にエクセルを使うことを想定して書かれているが、これをR（なりPythonなり）で書き直してみると、非常に勉強になると思う。

上の図に用いたコードを貼っておきます。

library(gganimate)
library(ggplot2)
library(tidyr)

x <- seq(0,1, by=0.01)
p <- 0.05
n <- length(x)
y <- rbinom(n,1,p)

plot(x,y,type = "h")

sim_pois1 <- function(t_n, p_n, w, seed){
  set.seed(seed)
  l <- 5 - 0.5*w
  r <- 5 + 0.5*w
  counts <- integer(t_n)
  p_list <- vector("list", t_n)
  for(i in 1:t_n){
    x <- runif(p_n, 0, 10)
    y <- runif(p_n, 0, 10)
    flag <- l < x & x < r & l < y & y < r
    p_list[[i]] <- data.frame(x=x,y=y,hit=flag,time=i)
    counts[i] <- sum(flag)
  }
  return(list(points=do.call("rbind",p_list),
              counts=counts))
}

sim_pois2 <- function(t_n, p, g_n, seed){
  set.seed(seed)
  counts <- integer(t_n)
  g_list <- vector("list", t_n)
  g <- expand.grid(x=seq(0,1,length.out=g_n),
                   y=seq(0,1,length.out=g_n))
  N <- nrow(g)
  for(i in 1:t_n){
    hit <- rbinom(N,1,p)
    g$hit <- hit
    g$time <- i
    g_list[[i]] <- g
    counts[i] <- sum(hit)
  }
  return(list(grid=do.call("rbind",g_list),
              counts=counts))
}


out <- sim_pois1(500,100,1,1)
out$points$facet_dummy <- "raw"
dfc <- do.call("rbind",lapply(1:500, function(i){
  tab <- table(out$counts[1:i])
  data.frame(time=i,
             num=as.integer(names(tab)),
             prob=as.integer(tab)/sum(tab),
             facet_dummy = "prob"
  )}))

dfseg <- data.frame(matrix(c(4.5,4.5,4.5,5.5,
                             4.5,4.5,5.5,4.5,
                             5.5,4.5,5.5,5.5,
                             4.5,5.5,5.5,5.5),4,4,byrow = TRUE),
                    facet_dummy = "raw")

rs <- range(out$counts)
dfprob <- data.frame(num=rs[1]:rs[2],
                      prob=dpois(rs[1]:rs[2],1),
                      facet_dummy = "prob")
ggplot()+
  geom_point(data=out$points, aes(x=x,y=y,colour=hit))+
  scale_color_viridis_d()+
  geom_col(data=dfc, aes(x=num, y=prob), width=0.7, position = "identity")+
  geom_point(data=dfprob, aes(x=num, y=prob), colour="grey")+
  geom_line(data=dfprob, aes(x=num, y=prob), colour="grey", linetype=2)+
  geom_segment(data=dfseg, aes(x=X1,y=X2,xend=X3,yend=X4), linetype=2)+
  facet_wrap(~facet_dummy, scales="free")+
  transition_time(time)+
  labs(title = 'Step: {frame_time}') +
  theme_bw(12)

anim_save("pois1.gif")

######

out2 <- sim_pois2(500,0.01,10,2)
out2$grid$facet_dummy <- "raw"

dfc2 <- do.call("rbind",lapply(1:500, function(i){
  tab <- table(out2$counts[1:i])
  data.frame(time=i,
             num=as.integer(names(tab)),
             prob=as.integer(tab)/sum(tab),
             facet_dummy = "prob"
  )}))



rs <- range(out2$counts)
dfprob2 <- data.frame(num=rs[1]:rs[2],
                     prob=dpois(rs[1]:rs[2],1),
                     facet_dummy = "prob")

ggplot()+
  geom_tile(data=out2$grid, aes(x=x,y=y,fill=hit))+
  scale_fill_viridis_c()+
  geom_col(data=dfc2,aes(x=num,y=prob), width=0.7, position = "identity")+
  geom_point(data=dfprob2, aes(x=num,y=prob), colour="grey")+
  geom_line(data=dfprob2, aes(x=num,y=prob), colour="grey", linetype=2)+
  facet_wrap(~facet_dummy, scales="free")+
  transition_time(time)+
  labs(title = 'Step: {frame_time}') +
  theme_bw(12)

anim_save("pois2.gif")

（このコードは見栄えを良くするためにやっていることが多く、無駄に難しくなってしまっているので、ポアソン³の想定読者のためにはもっとラフに書いたほうがよかったかもしれない。）

余談：最近ははてなブログよりZennを使っていて、Zennのほうが数式が書きやすくていい（マークダウン記法とTeX記法の衝突が少ない）のだけど、サイズの大きいGIFアニメは貼れないことがあるみたいです。

変な形の事後分布のアニメーション

ハイパーパラメータによって事後分布の形が大きく変わる様子です.

変な形の尤度関数をプロットする - ジョンとヨーコのイマジン日記 の続き的な投稿.

例1: 混合正規分布

モデルとして次の分布を考える.

ここで は平均0, 分散1 の正規分布の密度関数, とする.

さらに事前分布として, パラメータ にベータ分布 , に適当な幅の一様分布（密度が定数）を設定する.

この手のモデルは2つの分布が混ざったような形になるので混合分布とよばれ, モデルに基づくクラスタリングなどの際によく使われる.

ただし未知パラメータ（データから推定したいパラメータ）を2個にしてるのは単に図を作りやすくする都合のほうが大きい.

サンプルの大きさを100として, このときの事後分布をプロットしてみる. ごとに, 最大値が1（黄色っぽい色）になるようスケーリングしてある.

が小さいときはサポートの端っこに密度が張り付いていたようになっていたのが大きくなると真ん中によってくる様子.

例2: 非線形回帰

モデルとして次の分布を考える.

は既知の数とする.

事前分布として, パラメータ に正規分布 を設定する.

例えば を時間として, なんらかの施策の前後でデータの変化を見たいときは（データを適当にスケーリングするとして）このようなモデルを考えるかもしれない.

ちなみに一番単純なニューラルネットワークではこのような非線形の変換をベクトルに対して繰り返して行い, 例えば のようなモデルを考える.

例1のときと同様, サンプルの大きさを100として, このときの事後分布をプロットしてみる.

R のコード

まとめて貼ります.

library(dplyr)
library(ggplot2)
library(gganimate)

logsumexp2 =function (logx1,logx2){
  logx1 + log1p(exp(logx2-logx1))
}

logsumexp <- function(x){
  mx <- max(x)
  mx + log(sum(exp(x-mx)))
}

llmixnorm <- function(par, y){
  a0 <- par[1]
  b0 <- par[2]
  theta <- par[3]
  sum(logsumexp2(log(1-a0)+dnorm(y,log = TRUE), log(a0)+dnorm(y,b0,log = TRUE))) + dbeta(a0, theta, theta, log = TRUE)
}

llcp <- function(par, y, x){
  a0 <- par[1]
  b0 <- par[2]
  theta <- par[3]
  sum(dnorm(y,a0*tanh(b0*x),log = TRUE)) +
    dnorm(a0, 0, theta, log = TRUE)+
    dnorm(b0, 0, theta, log = TRUE)
}

N <- 100L
set.seed(1)
y11 <- c(rnorm(N/2),rnorm(N/2,0.5))

ggplot()+
  geom_histogram(data=NULL, aes(x=y11), fill="grey70", bins=25)+
  theme_bw(14)+labs(x="y", y="count")

parms <- expand.grid(a=seq(0.01,0.99,length.out = 60),
                     b=seq(0,1,length.out = 60),
                     theta=seq(0.5,2,by=0.1))

l11 <- apply(parms, 1, llmixnorm, y=y11)
dfL11 <- data.frame(parms, value = exp(l11)) %>% 
  group_by(theta) %>% 
  mutate(value=value/max(value)) %>% 
  ungroup()

ggplot(dfL11,aes(x=b,y=a))+
  geom_point(aes(colour=value), pch=15, size=3.5)+
  scale_colour_continuous(type = "viridis")+
  labs(title = 'theta: {round(frame_time,1)}') +
  transition_time(theta)+
  theme_bw(14)+theme(legend.position = "none")

anim_save("post_mixnorm.gif")

####
#tanh

parms <- expand.grid(a=seq(-2,2,length.out = 60),
                     b=seq(-10,10,length.out = 60),
                     theta=seq(1.5,0.1,by=-0.1))

x <- seq(-1,1,length.out=N)
set.seed(1);y11 <- c(rnorm(N/2,-.5),rnorm(N/2,0.5))

ggplot(data = NULL)+
  geom_point(aes(x=x, y=y11))+
  theme_bw(14)

l11 <- apply(parms, 1, llcp, y=y11,x=x)
dfL11 <- data.frame(parms, value = exp(l11)) %>% 
  group_by(theta) %>% 
  mutate(value=value/max(value)) %>% 
  ungroup()


ggplot(dfL11,aes(x=b,y=a))+
  geom_point(aes(colour=value),pch=15, size=3.5)+
  scale_colour_continuous(type = "viridis")+
  labs(title = 'theta: {round(frame_time,2)}') +
  transition_time(theta)+
  theme_bw(14)+theme(legend.position = "none")

anim_save("post_tanh.gif")