r - R中的中心极限定理

Question

我想模拟中心极限定理以证明它，但我不确定如何在 R 中进行。我想创建 10,000 个样本，样本大小为 n(可以是数字或参数)，从我将选择的分布中(均匀分布、指数分布等)。然后我想在一个图中绘制(使用 par 和 mfrow 命令)原始分布(直方图)，所有样本均值的分布，均值的 Q-Q 图，以及在第 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 指出给定 i.i.d.来自具有均值和方差的分布的样本，样本均值(作为随机变量)的分布随着样本数 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))
}

为了生成图表，我们使用 ggplot2 和来自 here 的 Aaron 的 qqplot.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) 原始生成分布的方差)。