ポアソン分布の2つの起源；島谷（2017）から - ジョンとヨーコのイマジン日記

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

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

作者:島谷健一郎
近代科学社

Amazon

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

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

もう一つ：

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

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

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

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

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

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

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

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

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

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

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

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

ポアソン³は主にエクセルを使うことを想定して書かれているが、これを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")