是的，这是正确的方法。 让我们看看它可以去哪里。

快速起步

您通常可以通过从内到外的工作来潜入一个完整的模拟。R命令_

x <- rcauchy(10, location=0, scale=1)

生成一个样本$10$并将其存储在变量中x。下一步是计算感兴趣的统计数据。这就是您将原始输出存储在变量中的原因：因此您可以多次使用它而无需重新生成它。

mean(x)
median(x)
sd(x)

等等等等

您可以返回第一个计算 (of x)，重复它，然后重新计算其平均值和标准差（或任何其他统计数据）。在某种程度上，这很有趣，也很有启发性。但是，对于模拟，您将需要重复这些过程数百（或更多）次。这是您希望将您的工作封装成可靠、可重用、灵活的组件的时候。

创建模拟

R（以及许多现代计算平台）通过一次生成所有模拟数据来发挥最佳作用。让我们通过编写一段模块化代码——一个函数——来处理样本生成和模拟迭代（在旧语言中需要双循环）来尊重这种偏好。它的投入显然至少包括：

• n, 样本量
• N，模拟大小（样本数）。让我们默认为$1$用于测试目的。
• ...，要传递给的任何其他参数rcauchy。

开始了：

sim <- function(n=10, N=1, ...) {
  x <- matrix(rcauchy(n*N, ...), nrow=n)             # Each sample is in a separate column
  stats <- apply(x, 2, function(y) c(mean(y), sd(y)))# Compute sample statistics
  rownames(stats) <- c("Mean", "SD")
  return(stats)
}

为输出的行命名的原因是为了帮助我们阅读它，就像在这个具有三个迭代的模拟中一样。首先出现命令，然后出现它的输出。

sim（N=3，位置=0，比例=1）

          [,1]      [,2]       [,3]
Mean -2.837735  6.471259 0.04831808
SD    4.445549 17.837725 7.00943078

已安排将所有均值放在第一行，将标准差放在第二行。

有了这个，让我们通过设置随机数种子、为独立随机样本生成一些统计数据并检查它们来运行可重现的模拟。让我们绘制直方图而不是使用汇总统计来描述两行输出（均值和 sds）。

set.seed(17)
x <- sim(n=10, N=500, location=0, scale=1)
hist(x["Mean", ])
hist(x["SD", ])

（老实说，我每次发布的结果都尽量使用相同的种子，这样人们就会知道我没有玩弄种子来调整结果，使它们看起来更像我的预期！对于私人探索我确实改变了种子，或者不设置它们，这样我每次都能看到新的结果。）

使用、重用和扩展模拟代码

现在你可以玩了。

这些直方图看起来很糟糕。是因为没有足够的模拟吗？重新运行最后四行代码但更改$N$至$5000$， 说。它没有帮助。改变location和scale。还在迷茫吗？也许我们应该回到更熟悉的情况。如何生成正态分布样本？为此，让我们扩展我们的主力函数sim 并让分布本身也成为一个参数！

sim <- function(n=10, N=1, f=rcauchy, ...) {
  # `f` is any function to generate random variables.  Its first argument must
  # be the number of values to output.
  x <- matrix(f(n*N, ...), nrow=n)                    # Samples are in columns
  stats <- apply(x, 2, function(y) c(mean(y), sd(y))) # Compute sample statistics
  rownames(stats) <- c("Mean", "SD")
  return(stats)
}

所做的唯一更改是包含f=rcauchy在参数列表中并将对的一个引用替换rcauchy为f。让我们用正态分布试试这个：

set.seed(17)
x <- sim(n=10, N=5000, f=rnorm, mean=0, sd=1)
par(mfrow=c(1,2)) # Shows two plots side-by-side
hist(x["Mean", ], main="Normal Means")
hist(x["SD", ], main="Normal SDs")

这似乎奏效了。现在，您可以通过编辑该x <- sim(...)行并重新运行它来改变分布、样本大小、模拟大小和分布参数。在一两分钟内，您应该可以很好地了解分布参数的变化如何改变模拟结果。通过更多的探索（考虑绘制均值和 SD 的序列），您应该能够明白为什么 Cauchy 样本分布看起来如此混乱。

为了研究柯西分布（和其他长尾分布），一个很好的修改sim是在其输出中包含样本中位数。

后续分析

最后，“均方误差”（MSE）呢？这些通常用于将估算器与其正在估算的值进行比较。我建议存储（在变量中）任何将被多次引用的值。因此，例如，您可以像这样研究平均值的均方误差：

location <- 0
x <- sim(n=10, N=5000, f=rcauchy, location=location, scale=1)
mean((x["Mean", ] - location)^2)                              # Mean squared error

如果要进行大量计算，请考虑将 MSE 的计算存储在函数中。

模拟输出的其他评估——绘图、描述性统计等——也很容易计算，因为相关模拟输出的整个数组仍然可以存储在x.

更进一步：来源和原则

有关 中工作模拟的更多示例，请在R我们的网站上搜索[R] Simulation。其中许多将很好地移植到其他通用计算平台，如 Python、Matlab 和Mathematica。 将它们移植到 Stata 或 SAS 等专门的统计平台可能需要做更多的工作，但同样的原则适用：

• 从内到外开发模拟代码（这通常与良好的软件开发实践相反！）。

• 封装有用的代码块，例如仿真代码、计算 MSE 的代码，甚至是绘制或汇总仿真输出的代码。

• 尊重计算环境的偏好和特性以获得最佳性能（但不要让它们支配你的注意力：你这样做是为了了解统计数据，而不是关于RPython 或其他什么）。

• 扩展封装代码的功能，而不是创建和修改原始代码的副本，以处理您正在做的事情的微小变化。

• 一次一小步地进行扩展，而不是一开始就尝试编写一个无所不能的功能。

• 使用良好的命名和注释约定来记录代码中没有立即清楚的任何内容。

• 可视化输出。

• 玩你的模拟，向他们学习。为此，使它们尽可能易于使用并以合理的速度运行。