引导程序,蒙特卡洛

机器算法验证 r 自习 引导程序 蒙特卡洛
2022-03-17 22:36:55

我已将以下问题设置为作业的一部分:

设计并实施一项模拟研究,以检查 bootstrap 的性能,以获得单变量数据样本平均值的 95% 置信区间。您的实现可以在 R 或 SAS 中。

您可能想要查看的性能方面是置信区间覆盖率(即置信区间包含真实均值的比例)和蒙特卡洛变异(即模拟之间的置信上限和下限变化多少)'

有谁知道如何处理这方面的蒙特卡洛变化?我似乎什至无法制定算法或任何东西。它与蒙特卡洛积分有关吗?谢谢!

4个回答

引导程序和蒙特卡洛程序之间的这种混淆不断重复出现,所以也许这是解决它的好地方。R代码示例也可能有助于完成作业。)

考虑以下引导程序的实现R

boot <- function(x, t) { # Exact bootstrap of procedure t on data x
    n <- length(x)       # Must lie between 2 and 7 inclusive.
    if (n > 7) {
        stop("Sample size exceeds 7; use an approximate method instead.")
    }
    p <- c(n, 1:(n-1))
    a <- rep(x, n^(n-1))
    dim(a) <- rep(n, n)
    y <- as.vector(a)
    while (n > 1) {
        n <- n-1
        a <- aperm(a, p)
        y <- cbind(as.vector(a), y)
    }
    apply(y, 1, t)
}

快速查看将确认这是一个确定性计算:没有生成或使用随机值。(我将把它的内部工作细节留给有兴趣的读者自己去了解。)

的参数boot是数组中的一批数值数据x和对函数的引用t(可以完全像 一样应用于数组x)以返回单个数值;换句话说,t是一个统计量它生成所有可能的带有替换的样本x并将其应用于t每个样本,从而为每个这样的样本生成一个数字:简而言之,这就是引导程序。返回值是一个数组,表示样本精确引导分布tx

作为一个小例子,让我们引导样本x=的平均值c(1,3)

> boot(c(1,3), mean)
> [1] 1 2 2 3

确实有四个可能的大小为的样本替换为将它们全部生成(按刚刚列出的顺序)并应用于它们中的每一个。在这种情况下,计算平均值,结果分别为,如输出所示。2(1,3)(1,1)(1,3)(3,1)(3,3)boottt1223

你从这里去哪里取决于你想如何使用引导程序。有关引导程序的完整信息包含在此输出数组中,因此显示它通常是个好主意。引导标准偏差的示例(1,3,3,4,7)

hist(boot(c(1,3,3,4,7), sd))

SD的直方图

现在我们准备谈谈蒙特卡罗模拟。假设,比如说,我们要通过使用自举分布的上 95% 百分位数, 个样本中引导 SD 的 95% 置信上限。这个过程有什么属性?找出答案的一种方法是假设样本是从均匀分布中随机获得的。(应用程序通常会指出一个合理的分布假设可能是什么;在这里,我随意选择了一个计算简单但分析上不容易处理的假设。)我们可以通过获取这样的样本并计算 UCL 来模拟发生的情况:5

> set.seed(17)
> quantile(boot(runif(5, min=0, max=10), sd), .95)[1]
     95% 
3.835870 

这个特定随机样本的结果是 3.83587。这是确定的:如果您boot再次使用相同的数据集进行调用,答案将完全相同。但是随着不同的随机样本,答案会如何变化呢?通过重复此过程几次并绘制结果的直方图来找出:

> boot.sd <- replicate(100, quantile(boot(runif(5, min=0, max=10), sd), .95)[1])
> hist(boot.sd)

模拟直方图

如果我们进行另一组模拟,随机抽取的结果会有所不同,创建一个(稍微)不同的直方图——但与这个没有太大不同。我们可以放心地使用它来了解 SD 的引导 UCL 是如何工作的。作为参考,请注意均匀分布的标准偏差(跨越范围从,如此处指定)等于正如人们希望任何 UCL 值得其盐一样,直方图中的大多数值(四分之三,或 0.75)都超过了这一点:01010/122.887

> length(boot.sd[boot.sd >= 10/sqrt(12)]) / length(boot.sd)
[1] 0.75

但这与我们指定(并且希望)的名义 95% 相去甚远!这是模拟的一个价值:它将我们的希望与实际发生的事情进行比较。(为什么会出现差异?我相信这是因为引导 SD 不适用于非常小的样本。)

审查

  • Bootstrap 统计在概念上与任何其他统计(如均值或标准差)相同;他们只是往往需要很长时间来计算。(请参阅代码中的警告消息boot!)

  • 蒙特卡罗模拟可用于研究引导统计如何因获取样本的随机性而变化。在这种模拟中观察到的变化是由于样本的变化,而不是自举的变化。

  • (此处未说明)因为 bootstrap 统计可能需要大量计算(显然,对于大小为 n 的样本最多 ^计算),因此可以方便地近似bootstrap 分布。这通常是通过创建一个“黑盒”程序从真正的引导分布中随机获取一个值并重复调用该程序来完成的。集体输出近似于精确分布。由于黑盒中的随机性,近似值可能会发生变化——但这种变化是近似过程的产物。它不是(从概念上)引导过程本身固有的。nnn

bootstrap 是一种 Monte Carlo 技术,因为它涉及某种类型的随机抽样。如果您在同一组数据上运行两次引导程序,您将得到不同的答案。您在引导程序中使用的样本越多,您获得的变化就越少。

覆盖率是指来自相同假设抽样分布的不同数据集的结果变化。

我也不确定“蒙特卡洛变体”本身的确切含义。当然,应该可以查看迭代之间有多少变化,例如上限(或下限)的值,例如(提示)。也许他们只希望您为蒙特卡洛而不是引导程序这样做?不过,这不是我对练习的要求。最好只问这个短语是什么意思。

我对这个家庭作业的理解是,它要求你对某种统计技术进行蒙特卡罗模拟。也就是说,您模拟一堆随机数据集,将此技术应用于这些数据集,并存储数字以便以后以方便的方式进行汇总(均值、模拟概率等)

现在,所讨论的技术是引导程序,其中涉及蒙特卡罗模拟(除非,正如 whuber 所证明的,您被要求执行精确的引导程序,这是极不可能的)。因此,作为模拟的结果,您可能会存储样本均值、样本标准差以及自举获得的均值的置信区间限制。

我认为为了使这项研究变得有趣,值得从偏态分布(指数)中抽取一个小样本()。这样,基于 CLT 的标准置信区间可能会表现不佳:例如,对于需要在任一端丢失 5% 的 90% 置信区间,您会看到类似于 10% 和 2% 的结果。bootstrap CI 可能会显示出更好的性能(例如,6% 和 4%)。n=10