R中的中心极限定理

Question

我希望模拟中心极限定理以证明它，但我不知道如何在 R 中做到这一点。我想从分布中创建 10,000 个样本大小为 n（可以是数字或参数）的样本我会选择（均匀、指数等...）。然后我想在一个图中（使用 par 和 mfrow 命令）绘制原始分布（直方图）、所有样本均值的分布、均值的 QQ 图，以及在第 4 个图中（有四个，2X2 )，我不确定要绘制什么。你能帮我开始用 R 编程吗？我想一旦我有了模拟数据我应该没问题。谢谢你。

我的初步尝试如下，它太简单了，我什至不确定是否正确。

r = 10000;
n = 20;

M = matrix(0,n,r);
Xbar = rep(0,r);

for (i in 1:r)
{
  M[,i] = runif(n,0,1);
}

for (i in 1:r)
{
  Xbar[i] = mean(M[,i]);
}

hist(Xbar);

Answer 1

CLT 指出，给定来自具有均值和方差的分布的 iid 样本，样本均值（作为随机变量）的分布会随着样本数量的n增加而收敛为高斯分布。在这里，我假设您要生成r包含n每个样本的样本集以创建r样本均值的样本。一些代码如下：

set.seed(123) ## set the seed for reproducibility
r <- 10000
n <- 200      ## I use 200 instead of 20 to enhance convergence to Gaussian

## this function computes the r samples of the sample mean from the 
## r*n original samples
sample.means <- function(samps, r, n) {
  rowMeans(matrix(samps,nrow=r,ncol=n))
}

为了生成图，我们使用这里的ggplot2Aaronqqplot.data函数。我们也用gridExtra在一帧中绘制多个图。

library(ggplot2)
library(gridExtra)
qqplot.data <- function (vec) {
  # following four lines from base R's qqline()
  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")  
}

generate.plots <- function(samps, samp.means) {
  p1 <- qplot(samps, geom="histogram", bins=30, main="Sample Histogram")
  p2 <- qplot(samp.means, geom="histogram", bins=30, main="Sample Mean Histogram")
  p3 <- qqplot.data(samp.means)
  grid.arrange(p1,p2,p3,ncol=2)
}

然后我们可以使用这些具有均匀分布的函数：

samps <- runif(r*n)  ## uniform distribution [0,1]
# compute sample means
samp.means <- sample.means(samps, r, n))
# generate plots
generate.plots(samps, samp.means)

我们得到：

或者，使用均值 = 3的泊松分布：

samps <- rpois(r*n,lambda=3)
# compute sample means
samp.means <- sample.means(samps, r, n))
# generate plots
generate.plots(samps, samp.means)

我们得到：

或者，具有平均值 = 1/1的指数分布：

samps <- rexp(r*n,rate=1)
# compute sample means
samp.means <- sample.means(samps, r, n))
# generate plots
generate.plots(samps, samp.means)

我们得到：

请注意，样本均值直方图的均值看起来都Gaussians与原始生成分布的均值非常相似，无论是均匀分布、泊松分布还是指数分布，如 CLT 预测的那样（其方差将为 1/ (n=200) 原始生成分布的方差）。