2-1. 緑本とggplot2: カウントデータとポワソン分布

0. はじめに

統計モデリングを勉強したくて、Twitterのフォロワーさん方絶賛の久保先生の通称 " みどりぼん " を読んでみました。もちろんそのままでもわかりやすいのですが、データから自分でグラフを作っていくと、より理解が深まると思いました。せっかくなので本編では使われていなかった{ggplot2}を使って、グラフを描くコードを書いていろいろと復習してみました。まずは第二章の確率分布と統計モデルの最尤推定から。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

図2.5 種子数の観測データとポアソン分布の対応

まずはこのグラフから。ヒストグラムはある植物の種子数の頻度を表現したもので、折れ線グラフはポアソン分布です。種子数はカウントデータ、つまりは非負の整数でなので、その確率分布はポアソン分布に近似しますよ、ということを示しています。

ポアソン分布は確率分布の1つで、そのパラメーターは " 分布の平均 " のみです。まずはこのグラフを描いてみましょう。

準備

# tidyverseをinstall
library(tidyverse)

# 種子数のデータをダウンロードする
# dataに格納されます
load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/distribution/data.RData"))

> data
[1] 2 2 4 6 4 5 2 3 1 2 0 4 3 3 3 3 4 2 7 2 4 3 3 3 4 3 7 5 3 1 7 6 4 6 5 2 4 7 2 2 6 2 4 5 4 5 1 3 2 3

# tbl_dfに変換して、列名をSeedにしておく
Data <-
  tbl_df(data) %>%
  rename(Seed=value)

# それぞれの種子数を持つ植物の個体数を出現頻度
# 確率分布に対応させるため
Data_Ratio <-
  Data %>%
  group_by(Seed) %>%
  summarise(Ratio = n()/nrow(Data))

> head(Data_Ratio)
# A tibble: 6 x 2
   Seed Ratio
  <dbl> <dbl>
1     0  0.02
2     1  0.06
3     2  0.22
4     3  0.24
5     4  0.20
6     5  0.10

ポアソン分布を生成する

ポアソン分布のパラメーターは平均値なので、与えられたデータの種子数の平均値を求めます

mean(data)

[1] 3.56

つまり、平均が3.56のポアソン分布がこの種子植物の種子数に近似しそうということになります。

# dataの平均3.56をとるポワソン分布においてSeedが0となる確率、1となる確率...と算出する
data.frame( Seed = 0:9 ) %>% # 種子数は0から9までと仮定する
  mutate( Probability = dpois(Seed, lambda=3.56)) -> Poisson

head(Poisson)

  Seed Probability
1    0  0.02843882
2    1  0.10124222
3    2  0.18021114
4    3  0.21385056
5    4  0.19032700
6    5  0.13551282

ggplot2で描画する

# Data_Ratio: 実際の観測値
# df: ポワソン分布から生成した確率分布
ggplot(Data_Ratio, aes(x=Seed,y=Ratio))+
  geom_bar(stat="identity",fill="white",colour="black")+ # 種子数のヒストグラムを描く
  geom_line(data=Poisson,aes(x=Seed,y=Probability))+ # 平均3.56のポワソン分布を重ねる
  geom_point(data=Poisson,aes(x=Seed,y=Probability),size=4,color="gray")+
  ylab("Ratio / Probability")
# ggsave("Figure/2-1.png", plot=last_plot(), dpi=300, width=6, height=4)

これで図2-5が描けると思います。

ポワソン分布とは何か？

[tex: p(y|\lambda)=\frac{\lambda^yexp(-\lambda)}{y!}]

ポワソン分布とは？

• p(y|λ): 平均がλで有る時にポアソン分布に従う確率変数がyとなる
• 平均λのみがポワソン分布のパラメーターである
• yは0以上の整数をとる
  • つまりカウントデータとなる
• すべてのyについての確率変数の和は1
• 確率分布の平均はλ ( λ >= 0 )
• 分散と平均が等しい

ポアソン分布に従うと考える時

• データに含まれている値が非負の整数である
  • カウントデータ
• yiに加減はあるだろうが上限がよくわからない
• 平均と分散がだいたい等しい

図2.6 さまざまな平均(λ)のポワソン分布

このグラフはλを変化させるとどう変わるか？ということを示したグラフです。それぞれ、λが3.5、7.7、15.1のポワソン分布を示しています。平均が大きくなれば大きくなるほど、グラフは右にずれて、ピークの値が小さくなっていくのがわかります。ポワソン分布は確率分布ですから、グラフの積分値は常に1となるので、横に広がれば、高さは低くなることがイメージできます。

data.frame(Seed=0:20) %>% # 種子数を0から20と仮定する
  mutate(Lambda3.5=dpois(Seed, lambda=3.5), # それぞれのλのポワソン分布を発生させる
         Lambda7.7=dpois(Seed, lambda=7.7),
         Lambda15.1=dpois(Seed, lambda=15.1)) %>%
  gather(Lambda,Probability,-Seed) -> df2

ggplot(df2, aes(x=Seed, y=Probability))+
  geom_line(aes(color=Lambda))+
  geom_point()
# ggsave("Figure/2-2.png", plot=last_plot(), dpi=300, width=6, height=4)

息切れしてきたのでこの辺で。また気分転換に投稿します。