可視化で理解する中心極限定理（ベルヌーイ分布編）

こういう話がある。

• 可視化で理解する中心極限定理 #rstatsj

非常に素晴らしいので、指数分布じゃなくて、ベルヌーイ分布版をアニメーションにしてみた。
下図は

• ベルヌーイ分布(コイン投げでいうところの表が出る確率p=0.2)からサンプリングしたデータのヒストグラム（灰色）
• 平均値の推定値（理論値、青い密度描画、中心極限定理から平均0.2, 標準偏差sqrt(0.2*0.8/データ数)の正規分布に従う）
• データから推定した平均値、およびその１σ信頼区間（赤色縦棒、1σ信頼区間は正規分布の1σにそのまま対応）

となっている。

データが溜まってくると平均値の推定値の精度があがる（赤線の間隔 or 正規分布の標準偏差が縮まる）ことが可視化されてわかりやすい。

コードは以下で、ほぼ@hoxo_m氏の上記の記事のパクリ。ggplot2の勉強になりました。

library(animation)
library(ggplot2)
saveGIF( {
  prob <- 0.2
  mean_p <- prob
  sd_p <- sqrt(prob*(1-prob))
  x <- c()
  for(i in 1:50) {
    x <- c(x, sample(c(1,0), 10, prob=c(prob, 1-prob), replace=TRUE))
    mean_x <- mean(x)
    sd_x <- sd(x)/sqrt(length(x))
    p <- ggplot(data.frame(x=x), aes(x=x)) + 
      geom_histogram(aes(y=..density..), binwidth=0.1, alpha=0.5) + 
      stat_function(fun=dnorm, geom="density", color="blue", fill="blue", alpha=0.5, arg=list(mean=mean_p, sd=sd_p/sqrt(length(x)))) + 
      geom_vline(xintercept=mean_x,      color="red", size=2) +
      geom_vline(xintercept=mean_x+sd_x, color="red", size=1) +
      geom_vline(xintercept=mean_x-sd_x, color="red", size=1) +
      scale_y_continuous(labels = function(x) sprintf("%0.2f", x)) + 
      xlim(0, 1.1) + xlab("ベルヌーイ分布の平均値の分布")      
    print(p)
  }
}, cmd.fun = system, interval = 0.4, ani.width=600, ani.height=400)

推定した平均値のバラつき（正規分布の１σ or データから計算した1σ信頼区間）も、理論・データからの推計値がほぼ一致。

> sd_x
[1] 0.01823362
> sd_p/sqrt(length(x))
[1] 0.01788854