今回は、Data Visualization Advent Calendar 2015への寄稿として書きます。qiita.com
データを特徴づける指標にはいろいろなものがあります。例えば、算術平均は分布の中心を表す代表的な指標です。これに加えて、データの形状やばらつきを把握すると、より詳しくデータの性質を理解できます。今回は、データの形状やばらつきに関連する話題として中心極限定理（Central Limit Theorem）を取り上げます。特に、中心極限定理の内容を可視化することで、理解しやすくなるかどうかを試みます。
中心極限定理
まず中心極限定理の説明をしておきます。
この定理は、母集団から抽出したサンプルの性質を決めています。母集団からサンプルを抽出する場合、（１）ある程度のサンプルサイズがあれば、その平均値（もしくは和）の分布が正規分布（ガウス分布）になると決まっています。さらに、この定理は、（２）「母集団の形状が正規分布でなくても」成り立ちます。
中心極限定理の可視化
上の文章は、いかにも統計学のテキストに載っている感じなのですが、率直にいって「わかりにくい」です。抽象的な定理がわかりにくい時には、具体例で考えれば良いのですが、ちょっと前までは、それもそんなに簡単ではありませんでした。なぜかというと、実際にデータを取るような場合、まず母集団の全体像が把握できないからです。さらに、サンプルの集団を繰り返し抽出するといっても、それに時間や労力のコストがかかることも実例で試す難しさを上げてしまいます。しかし、最近は、個人のPCでもシミュレーションを簡単に作ることができるようになっています。そこで、中心極限定理に関する（１）と（２）をシミュレーションで確かめてみましょう。シミュレーションの構築にはR言語を使用します。
正規分布をとる母集団の場合
まず、正規分布を取る母集団を仮定します。わかりやすい例として平均０、標準偏差１の正規分布にします。そのような母集団からランダムに１００個の値を取り出すには、以下のように書きます。
sample <- rnorm(100, mean = 0, sd = 1)
さらにこのサンプルの平均値は以下のように計算します。
sampleMean <- mean(sample)
これで１回のサンプルの平均値が求まります。ただし、これはランダムに取り出したサンプルなので、取り出すたびに計算結果が変わります。そこで、この計算を１００００回繰り返して、サンプルの平均値の分布（サンプリング分布）を作ってみましょう。
#10000回分の平均値を保存するためのオブジェクト res <- rep(NA,10000) for (i in 1:10000){ sample <- rnorm(100, mean = 0, sd = 1) sampleMean <- mean(sample) res[i] <- sampleMean } #サンプルの平均値の分布をプロット hist(res)
左右対称の形が出てきました。この分布は、平均が母平均の期待値に一致します。また、標準偏差は、母集団の標準偏差をサンプル数の平方根でわったもの（いわゆる標準誤差）になっています*1。この例の場合は、平均０、標準偏差0.1となります。そのような正規分布曲線を実際に重ね合わせたものは以下のようなグラフになります。
hist(res) #次のグラフを重ね合わせるため設定 par(new = TRUE) #平均０、標準偏差０．１の正規分布曲線 curve(dnorm(x,mean =0, sd = 0.1),from = min(res), to = max(res), axes = FALSE, xlab = "",ylab = "")
正規分布曲線がシミュレーションで作成したヒストグラムの上をきれいに通り抜けています。
これで、サンプルの平均値が正規分布していることが確認できました。
正規分布以外をとる母集団の場合
上記で、もともと正規分布を取る母集団から抽出したサンプルの平均値が正規分布をするというのは、なんとなくそういうものかなと思えます。
では、正規分布以外のデータでも確認してみましょう。今回は、サイコロを振る場合を例として考えます。サイコロは１から６の目が等頻度で出てくるので一様分布を仮定できます。
R言語ではsample関数を利用すれば簡単に一様分布を持つ母集団をつくれます。
dice <- sample(1:6,100000,replace = TRUE) hist(dice, breaks = c(1:7), right = FALSE)
では１０回サイコロを振った平均値の分布を計算するシミュレーションを１００００回行って、分布を調べてみましょう。
#10000回分の平均値を保存するためのオブジェクト res <- rep(NA,10000) for (i in 1:10000){ dice <- sample(1:6,10,replace = TRUE) diceMean <- mean(dice) res[i] <- diceMean } hist(res)
先ほどの例と同じように左右対称の分布が出てきました。この分布の分散は一様分布の標準偏差をサンプルサイズの平方根で割ったものです。一様分布の分散は(最高値ー最低値)^2/12で計算できます。今回のサンプルサイズは１００なので、これらをまとめると標準誤差はおよそ0.144となります。平均3.5, 標準偏差0.144の正規分布曲線を書き足してみると以下の図のようになります。
まとめ：あらたな統計学の学び方としての可視化・シミュレーション
これまで、統計学を学ぶ時には、「定義を言葉で学ぶ」、「数式やその展開、証明から学ぶ」という二つの選択肢がありました。今日では、これに「データでシミュレーションをしながら学ぶ」という第３の選択肢が加わったように思います。統計学のテキストで述べられている定理、例題をどんどんシミュレートして、統計に関する感覚を養うといいと思います。
おまけ：「インタラクティブな」可視化で理解する中心極限定理
RのShinyパッケージを利用すると、平均やサンプルするなどを可変にしたシミュレーションをブラウザ上に表示できます*2。今回の中心極限定理を理解するためのシミュレーションのコードは以下の通りです。母集団として正規分布か一様分布を選べて、抽出するサンプルサイズを変えられます。特にサンプルサイズによって、サンプル平均のばらつきが変化することを理解できます。
shinyパッケージではui.Rとsever.Rという二つのファイルを使用します。u.Rとserver.Rという名前をつけたファイルに以下のそれぞれをコピーして、同一のフォルダ内に保存してください。またshinyパッケージをインストールしていない場合は、あらかじめインストールしておいてください。
ui.R
library(shiny) shinyUI(fluidPage( titlePanel("可視化で理解する中心極限定理"), sidebarLayout( sidebarPanel( radioButtons("dist1", "Distribution type:", c("Normal" = "norm", "Uniform" = "unif")), br(), sliderInput("n", "Number of samples:", value = 50, min = 6, max = 200) ), # Show a tabset that includes a plot, summary, and table view # of the generated distribution mainPanel( tabsetPanel(type = "tabs", tabPanel("Plot", plotOutput("plot")), tabPanel("AveragePlot",plotOutput("aveplot")) ) ) ) ))
server.R
library(shiny) shinyServer(function(input, output) { data1 <- reactive({ dist1 <- switch(input$dist1, norm = rnorm, unif = runif, rnorm) if(input$dist1 == "unif"){ (dist1(10000) - 0.5) * 6 } else { dist1(10000) } }) output$plot <- renderPlot({ dist1 <- input$dist1 n <- 10000 hist(data1(), main=paste('r', dist1, '(', n, ')', sep='')) }) output$aveplot <- renderPlot({ n <- input$n data2 <- rep(0,10000) for (i in 1:length(data2)){ data2[i]<- mean(sample(data1(),input$n)) } hist(data2, main=paste('average of ',n, ' samples ', sep='')) }) })
