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

想像してください。「あなたはぼくをプラグマティストだと言うかもしれない」と歌う、逆イマジンです。

R: 散布図の点が多すぎると感じるときに試すこと

散布図の点の数が多すぎると次のような問題を感じることがある。

  • 点が重なりまくる
  • 描画が遅い
library(ggplot2) #ggplot2 の diamonds データを例に使う
plot(price~carat, data=diamonds)

これらへの簡単な対処法として、pch="." を指定して点を小さくするというのがある。

plot(price~carat, data=diamonds, pch=".")

点の重なりが多少改善され、描画も少し速くなる(気がする)。

プロットに多少時間がかかるのはしようがないとして、マーカーの透明度を上げるという手もある。

#rgb() の4つめの引数が不透明度
plot(price~carat, data=diamonds, col=rgb(0,0,0,0.02))

点が重なったところは色が濃くなるので頻度の情報も多少読み取りやすくなる。

この2つの合わせ技が次の図:

plot(price~carat, data=diamonds, col=rgb(0,0,0,0.1), pch=".")

ggplot2 では、透明度を上げるには次のようにする。

theme_set(theme_bw(14)) #デフォルトのグレー背景が好みでないので
ggplot(diamonds, aes(carat, price, colour=cut))+
  geom_point(alpha=0.1)

ついでに色分けで層別してあるが、凡例の色も薄くなってしまって見えづらいときがある。

そんな場合のために guides 関数で凡例の alpha を上書きすることができる。

ggplot(diamonds, aes(carat, price, colour=cut))+
  geom_point(alpha=0.1)+
  guides(colour = guide_legend(override.aes = list(alpha = 1)))

このやり方は からっぽのしょこid:anemptyarchive さんに教えてもらった。ありがとうございます!

さらに点が多いときなどは hexbin を使うことも多い。

ggplot(diamonds, aes(carat, price))+
  geom_hex()

ビンごとに集計して度数を色の濃さ(など)で表すもの。

散布図の描画が遅いのは点の座標をデータの数だけポチポチしているからなので(だと思う)、集計してしまうと描画も速くなるし頻度も読みとりやすい。

逆にデメリットとしてはビンの幅によってデータを均しすぎたりもできるという可能性はあるかもしれない。

飛び抜けて値の大きいビンがあるとき(0が多いデータなど)は色を対数スケールにしたりすることもある。

ggplot(diamonds, aes(carat, price))+
  geom_hex(aes(fill=after_stat(log2(count))))

色分けと hexbin の合わせ技が次の図:

ggplot(diamonds, aes(carat, price, colour=cut))+
  geom_hex(aes(alpha=after_stat(log2(count))), fill="darkgrey")

でも素直に facet_ で図を並べたほうが見やすいかも:

ggplot(diamonds, aes(carat, price))+
  geom_hex(aes(alpha=after_stat(log2(count))), fill="darkgrey")+
  facet_wrap(~cut)

以上です。

ポアソン分布の2つの起源;島谷(2017)から

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

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

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

もう一つ:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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