r - R中的中心极限定理

我想模拟中心极限定理以进行演示，但不确定如何在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);

最佳答案

CLT指出了给定的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）原始生成分布的方差）。