중심극한정리

# ref : https://stackoverflow.com/questions/40307510/central-limit-theorem-in-r
# ref : https://ko.wikipedia.org/wiki/%EC%A4%91%EC%8B%AC%EA%B7%B9%ED%95%9C%EC%A0%95%EB%A6%AC
#
#  확률분포를 가진 독립 확률 변수 n개의 [평균의 분포]는 n이 적당히 크다면 정규분포에 가까워진다는 정리

library(ggplot2)
library(gridExtra)

Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

The following object is masked from 'package:randomForest':

    combine

# 모집단 : N(0, 1) , 10000000 elements
population = rnorm(n = 10000000, m = 0, sd = 1)
plot(density(population), main = paste("mean = ", mean(population), "/ sd = ", var(population)))

# 평균들의 분포 
sample.means <- function(population, sample_size, repeat_count) {
  rowMeans(matrix(population, nrow = repeat_count , ncol = sample_size))
}
qqplot.data <- function (vec) {
  y <- quantile(vec[!is.na(vec)], c(0.25, 0.75))
  x <- qnorm(c(0.25, 0.75))
  slope <- diff(y)/diff(x)
  int <- y[1L] - slope * x[1L]
  d <- data.frame(resids = vec)
  ggplot(d, aes(sample = resids)) + stat_qq() + geom_abline(slope = slope, intercept = int, colour="red") + ggtitle("Q-Q plot")  
}
eval <- function(population, sample_size, repeat_count) {
  sample.means.small <- sample.means(population, sample_size, repeat_count)
  simul_variance <- round(var(sample.means.small),5)
  expected_variance <- round(var(population) / sample_size,5)
  diff <- round(abs(expected_variance - simul_variance),5)
  
  p1 <- qplot(sample.means.small, 
        geom="histogram", 
        bins=30, 
        main=paste("sample_n : ", sample_size ," variance : ", simul_variance, " vs ", expected_variance, "(", diff , ")  mean : ", round(mean(sample.means.small),5)  ))
  p2 <- qqplot.data(sample.means.small)
  
  grid.arrange(p1,p2,nrow=2)
}
eval(population, 10, 1000)

eval(population, 50, 1000)

eval(population, 200, 1000)

# 모집단 : chi squared distribution , 10000000 elements
population = rchisq(n = 10000000, df = 0, ncp = 2.)
plot(density(population), main = paste("mean = ", mean(population), "/ sd = ", var(population)))

eval(population, 10, 1000)

eval(population, 50, 1000)

eval(population, 200, 1000)