练习中使用循环的说明是不好的建议。该rbinom函数已经能够模拟值向量，因此不需要循环。这里要做的最简单的事情是创建一个$r \times N$模拟伯努利随机变量矩阵$N=100$使您有足够的样本量来满足四个指定值的要求$n$. 假设您不介意嵌套模拟（这不是问题），那么您将拥有构建图形所需的所有模拟值。这是一些简单的代码来生成模拟伯努利随机变量的“可重现”矩阵。（请注意，您也可以使用该sample.int函数有效地模拟伯努利随机变量。）它不会为您提供该书中图表中使用的确切结果，但它仍然会为您提供可重现的模拟值。

#Set parameters
N    <- 100
r    <- 10000
PROB <- 0.78

#Create matrix for simulated values
set.seed(1)
SIMULATIONS <- matrix(rbinom(r*N, size = 1, prob = PROB), nrow = r, ncol = N)
colnames(SIMULATIONS) <- sprintf('Sample[%s]', 1:N)

一旦你模拟了矩阵SIMULATIONS，你就有了$r$模拟值的行$N$样本值。rowMeans您可以通过对矩阵的相关子集使用函数来获得样本均值的相关模拟。然后，您可以使用适当的绘图函数来构建所需的绘图。这给出了与您显示的图表相似的结果。

#Create matrix of standardised sample means
STD.MEANS <- matrix(0, nrow = R, ncol = 4)
colnames(MEANS) <- c('n[2]', 'n[5]', 'n[25]', 'n[100]')
STD.MEANS[, 1] <-   sqrt(2)*(rowMeans(SIMULATIONS[, 1:2])   - PROB)/sqrt(PROB*(1-PROB))
STD.MEANS[, 2] <-   sqrt(5)*(rowMeans(SIMULATIONS[, 1:5])   - PROB)/sqrt(PROB*(1-PROB))
STD.MEANS[, 3] <-  sqrt(25)*(rowMeans(SIMULATIONS[, 1:25])  - PROB)/sqrt(PROB*(1-PROB))
STD.MEANS[, 4] <- sqrt(100)*(rowMeans(SIMULATIONS[, 1:100]) - PROB)/sqrt(PROB*(1-PROB))

#Plot the histograms
par(mfrow = c(2,2))
hist(STD.MEANS[, 1], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 2)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)
hist(STD.MEANS[, 2], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 5)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)
hist(STD.MEANS[, 3], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 25)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)
hist(STD.MEANS[, 4], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 100)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)

首先，我同意@Ben(+1) 关于尽可能避免显式循环的声明。我使用for了循环，因为它们似乎是您的锻炼所必需的。

标准化是在for循环外完成的，使用 10,000 个平均值的平均值和标准差a。

这是R中的模拟案例$n = 25.$

set.seed(121)
n = 25;  p = 0.78
r = 10^4;  a = numeric(r)
for(i in 1:r) {
  a[i] = mean(rbinom(n, 1, .78))  
  }
mean(a);  sd(a)
z = (a-mean(a))/sd(a)
cp = seq(-5.75, 5.75, length=13)
hdr = "n=25: Standardized Value of Sample Average"
hist(z, prob=T, br=cp, ylim=c(0,.4), col="skyblue2", main=hdr)
 curve(dnorm(x), add=T, lwd=2)

关于制作直方图的注意事项：我使用了 R 基础的标准图形来制作直方图并叠加标准正态密度曲线。虽然$r = 10,000$手段已经生成，所以a 有$r$值，在我的模拟中 --fifteen 中没有很多唯一值（有些相对罕见）。a如果并列值在直方图箱中的比例不均等，您会得到一些看起来很奇怪的直方图，这会使 的值z看起来很不正常。通过选择 13 个箱子，我得到了一个不错的情节。（粗略地说，每个 bin 有两个z- 值，末端有一些空 bin。）

length(unique(a))
[1] 15
table(a)
a
0.44 0.48 0.52 0.56  0.6 0.64 0.68 0.72 0.76  0.8 0.84 0.88 0.92 0.96    1 
   1   12   22   86  212  459  847 1381 1795 1914 1591 1037  495  123   25