中心极限定理模拟

一、理论阐述

中心极限定理

设随机变量$X_1,X_2,...,X_n,...$独立同分布，并存在有限的期望与方差$E({X_i})=\mu,D({X_i})=\sigma^2(i=1，2，...)$，则对任意$x$，分布函数

$$F_n(x)=P\left\{\frac{\sum_{i=1}^n X_i-n \mu}{\sigma \sqrt{n}} \leq x\right\}$$

满足

$$\lim _{n \rightarrow \infty} F_n(x)=\lim _{n \rightarrow \infty} P\left\{\frac{\sum_{i=1}^n X_i-n \mu}{\sqrt{n} \sigma} \leq x\right\}=\frac{1}{\sqrt{2 \pi}} \int_{-\infty}^x e^{-\frac{t^2}{2}} d t$$

中心极限定理推论

从上述定理可得，对随机变量$X_1,X_2,...,X_n,...$，存在$E({X_i})=\mu,D({X_i})=\sigma^2(i=1，2，...)$，当$n$足够大时，满足

$$\frac{\sum\limits_{k=1}^nX_k-n\mu}{\sqrt{n}\sigma}\sim N(0,1)$$

由于

$$\frac{\sum_{\mathrm{k}=1}^{\mathrm{n}} \mathrm{X}_{\mathrm{k}}-\mathrm{n} \mu}{\sqrt{\mathrm{n}} \sigma}=\frac{\frac{1}{\mathrm{n}} \sum_{\mathrm{k}=1}^{\mathrm{n}} \mathrm{X}_{\mathrm{k}}-\mu}{\sigma / \sqrt{\mathrm{n}}}=\frac{\overline{\mathrm{X}}-\mu}{\sigma / \sqrt{\mathrm{n}}}$$

故有

$$\frac{\overline{\mathrm{X}}-\mu}{\sigma / \sqrt{\mathrm{n}}} \sim N(0,1),即\overline{\mathrm{X}} \sim N(\mu,\sigma^2/n)$$

二、算法

根据上述理论，我们可以得到一个通用的中心极限定理模拟算法

1. 生成$N(N \gg n)$个随机数，将其作为后续算法抽样的集合，记为$A$

2. 从集合$A$中抽取$n$个样本，并计算其均值，记为$\overline{\mathrm{X}}_1$

3. 将上述步骤重复$k$次，得到随机变量$\overline{\mathrm{X}}_1,\overline{\mathrm{X}}_2,...,\overline{\mathrm{X}}_k$

4. 对上述随机变量进行检验，判断其是否符合$N(\mu,\sigma^2/n)$

三、R程序

rm(list=ls())
library(moments)
set.seed(0721)  # 设置随机数种子
# 定义函数，用于模拟中心极限定理
Central_Limit_Theorem <- function(pop, m, n = c(1, 5, 20, 50)){
  # pop：随机数  m:抽样次数
  mean_pop <- mean(pop) # 总体均值
  var_pop <- var(pop) # 总体方差
  # 随机抽取不同大小的样本，并计算样本均值的分布
  x_bar <- numeric(m) # 初始化一个向量存储样本均值
  par(mfrow = c(2, 4)) # 设置绘图区域为2*4
  # 初始化，用于存储描述性统计的矩阵
  result <- matrix(0, nrow = 4, ncol = 7)
  k <- 0
  for (i in n) {
    k <- k + 1
    for (j in 1:m) {
      x_bar[j] <- mean(sample(pop, size = i)) # 随机抽取i个样本，并计算样本均值
    }
    # 描述性统计
    result[k, ] <- c(mean_pop, mean(x_bar), var_pop, var(x_bar), var_pop/i, skewness(x_bar), kurtosis(x_bar))
    # K-S检验
    ks <- ks.test(x_bar, "pnorm", mean = mean (x_bar), sd = sd (x_bar))
    hist(x_bar, freq = FALSE, main = paste("n =", i), xlab = "Sample mean") # 绘制直方图
    curve(dnorm(x, mean = mean(x_bar), sd = sd(x_bar)), add = TRUE, col = "red")
    # Q-Q图
    qqnorm (x_bar, main = paste("p.value =", round(ks$p.value, 3))) # 绘制Q-Q图
    qqline (x_bar) # 添加参考线
  }
  # 创建一个数据框
  M <- data.frame(mean_X = c(result[, 1]), x_mean = c(result[, 2]), var_X = c(result[, 3]), x_var = c(result[, 4]), var_X_n = c(result[, 5]), skewness = result[, 6], kurtosis = result[, 7])
  return(M)
}

# 标准正态分布
pop <- rnorm(10000, mean = 0, sd = 1) # 总体大小为10000，均值为0，标准差为1
M <- Central_Limit_Theorem(pop, 1000, c(1, 5, 10, 30))

# 指数分布，rate为1
pop <- rexp(10000, rate = 1)
M <- Central_Limit_Theorem(pop, 1000)

# 均匀分布，假设最小值为0，最大值为1
pop <- runif(10000, min = 0, max = 1)
M <- Central_Limit_Theorem(pop, 1000)

# 泊松分布，假设lambda为1
pop <- rpois(10000, lambda = 1)
M <- Central_Limit_Theorem(pop, 1000, c(1, 10, 30, 70))

# 二项分布，假设试验次数为10，成功概率为0.5
pop <- rbinom(10000, size = 10, prob = 0.5)
M <- Central_Limit_Theorem(pop, 1000)

# 两点分布，假设两个可能的值为-1和1，概率相等
pop <- sample(c(-1, 1), size = 10000, replace = TRUE)
M <- Central_Limit_Theorem(pop, 1000, c(1, 30, 100, 200))

# write.csv(M, "ciallo.csv", row.names = TRUE)  # 将数据存入Excel

四、结果分析

分析方法

对于各个分布，每次抽样取$n$个样本用于计算均值，共取1000次，对于每个分布分别取4个不同的$n$值进行模拟，并各自进行检验，检验方法如下：

1. 直方图与核密度估计曲线
2. Q-Q图
3. K-S检验，检验其是否符合理论上的正态分布（p值绘制于Q-Q图上方）
4. 描述性统计，包括分布均值方差，理论方差，样本的均值，方差

正态分布

直方图，Q-Q图与K-S检验结果

描述性统计

n	原分布均值	样本均值	原分布方差	样本方差	理论样本方差
1	-0.0186	-0.0293	0.9902	0.9663	0.9902
5	-0.0186	-0.0092	0.9902	0.2081	0.198
10	-0.0186	-0.0313	0.9902	0.1015	0.099
30	-0.0186	-0.0125	0.9902	0.0343	0.033

指数分布

直方图，Q-Q图与K-S检验结果

描述性统计

n	原分布均值	样本均值	原分布方差	样本方差	理论样本方差
1	0.9998	0.9586	0.9822	0.9291	0.9822
5	0.9998	1.0008	0.9822	0.1954	0.1964
20	0.9998	1.0013	0.9822	0.0503	0.0491
50	0.9998	1.0002	0.9822	0.0195	0.0196

均匀分布

直方图，Q-Q图与K-S检验结果

描述性统计

n	原分布均值	样本均值	原分布方差	样本方差	理论样本方差
1	0.5001	0.5018	0.0842	0.0831	0.0842
5	0.5001	0.5026	0.0842	0.0165	0.0168
20	0.5001	0.5004	0.0842	0.0041	0.0042
50	0.5001	0.5007	0.0842	0.0016	0.0017

泊松分布

直方图，Q-Q图与K-S检验结果

描述性统计

n	原分布均值	样本均值	原分布方差	样本方差	理论样本方差
1	1.0137	0.992	1.0192	0.9609	1.0192
10	1.0137	1.0205	1.0192	0.1041	0.1019
30	1.0137	1.0062	1.0192	0.0334	0.034
70	1.0137	1.0077	1.0192	0.0143	0.0146

二项分布

直方图，Q-Q图与K-S检验结果

描述性统计

n	原分布均值	样本均值	原分布方差	样本方差	理论样本方差
1	4.9825	4.981	2.4702	2.379	2.4702
5	4.9825	4.9592	2.4702	0.4893	0.494
20	4.9825	4.977	2.4702	0.1198	0.1235
50	4.9825	4.9865	2.4702	0.0519	0.0494

两点分布

直方图，Q-Q图与K-S检验结果

描述性统计

n	原分布均值	样本均值	原分布方差	样本方差	理论样本方差
1	0.0146	-0.006	0.9999	1.001	0.9999
30	0.0146	0.0073	0.9999	0.0311	0.0333
100	0.0146	0.0119	0.9999	0.0104	0.01
200	0.0146	0.0132	0.9999	0.0054	0.005

总结

从上述结果来看，随着$n$的增大，各个样本会越来越趋近于正态分布，从直方图与核密度估计曲线可以看到直方图与曲线越来越接近，而且Q-Q图上的点也逐渐分布于参考线上，且在$n$值足够大时，所有的分布均可以通过K-S检验。而从描述性统计中可以看出，随着$n$的增大，样本的均值越来越接近分布均值，方差也越来越接近理论方差。