大数の法則と中心極限定理を実感する

一様分布からのランダムサンプリング

一様分布\(U(0,1)\)から，サンプルサイズ\(n=1,10,100,1000\)のサンプルをとって標本平均を計算．これを，それぞれ10000回繰り返す．赤線は中心極限定理による標本分布近似の理論的予想．

repdata1<-replicate(10000,mean(runif(1,0,1)))
repdata2<-replicate(10000,mean(runif(10,0,1)))
repdata3<-replicate(10000,mean(runif(100,0,1)))
repdata4<-replicate(10000,mean(runif(1000,0,1)))
hist(repdata1, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")

hist(repdata2, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")
curve(dnorm(x, mean=1/2, sd= sqrt(1/12)/sqrt(10)), type="l", col = "red", add=T)

hist(repdata3, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")
curve(dnorm(x, mean=1/2, sd= sqrt(1/12)/sqrt(100)), type="l", col = "red", add=T)

hist(repdata4, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")
curve(dnorm(x, mean=1/2, sd= sqrt(1/12)/sqrt(1000)), type="l", col = "red", add=T)

ベータ分布\(Beta(0.5,0.5)\)からのランダムサンプリング

repdata5<-replicate(10000,mean(rbeta(1,0.5,0.5)))
repdata6<-replicate(10000,mean(rbeta(10,0.5,0.5)))
repdata7<-replicate(10000,mean(rbeta(100,0.5,0.5)))
repdata8<-replicate(10000,mean(rbeta(1000,0.5,0.5)))
hist(repdata5, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")

hist(repdata6, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")
curve(dnorm(x, mean=1/2, sd= sqrt(1/8)/sqrt(10)), type="l", col = "red", add=T)

hist(repdata7, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")
curve(dnorm(x, mean=1/2, sd= sqrt(1/8)/sqrt(100)), type="l", col = "red", add=T)

hist(repdata8, prob=TRUE, breaks = seq(0,1,0.005), xlim=c(0,1), col="skyblue")
curve(dnorm(x, mean=1/2, sd= sqrt(1/8)/sqrt(1000)), type="l", col = "red", add=T)

おまけggplotで書く

要tidyverse

library(tidyverse)

## -- Attaching packages --------------------------------------------------------- tidyverse 1.2.1 --

## √ ggplot2 3.2.0     √ purrr   0.3.2
## √ tibble  2.1.2     √ dplyr   0.8.1
## √ tidyr   0.8.3     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.4.0

## -- Conflicts ------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

tibble(n1=repdata1,n10=repdata2,
           n100=repdata3,n1000=repdata4) %>% 
  gather(key=n,value=value,n1,n10,n100,n1000) %>% 
  ggplot() +
  geom_histogram(aes(x = value, y =..density..),bins=200,fill="skyblue", color= "black") +
  # geom_density(adjust = 4, colour='black', size=1.25) +
  facet_grid(n~.)