这是蒙特卡洛模拟吗？

机器算法验证 r 假设检验 蒙特卡洛
2022-04-10 16:24:46

所以，让我们比较两个正态分布

Do this x times: 

runs <- 100000
a.samples <- rnorm(runs, mean = 5) 
b.samples <- rbeta(runs, mean = 0) 
mc.p.value <- sum(a.samples > b.samples)/runs

低于我们的 alpha (0.05) 的 mc.p.values 除以 x 将给出 type1 错误率。我们的 H0 是 a.samples >= b.samples。（灵感来自https://www.countbayesie.com/blog/2015/3/3/6-amazing-trick-with-monte-carlo-simulations）

但是，我认为蒙特卡罗模拟必须遵循以下步骤：

算法：

  1. 为数据设置一些分布，f() 或 f(θ)，以及一些 H0
  2. 多次重复以下两个步骤： (a) 根据 H0 模拟一个数据集 (b) 使用模拟数据计算 T(x)
  3. 添加从样本数据评估的 T(X)
  4. 订购所有的 T(x)s
  5. p 值是 T(x) 与样本数据中的值相比极端或更极端的比例

因此，第一个代码片段不是真正的蒙特卡罗模拟吗？并且 p 值是否有效，因为如果您将其绘制成图表，您不会得到预期的 5% type1 错误率，这是人们可能期望的统计测试。

1个回答

这是一个蒙特卡洛模拟，虽然我怀疑它正在做你想做的事，而且真的没有实际用处。您要比较的是 10000 个单样本研究，并确定这些单独观察值中有多少平均更高。因此，最好将您的代码概念化为以下效率较低的代码：

runs <- 10000
result <- logical(runs)

for(i in 1:runs){
  a <- rnorm(1, mean = 0) 
  b <- rnorm(1, mean = 0)
  result[i] <- a > b
} 
mc.p.value <- mean(result)

上面的代码，当分布相等时，应该表明a高于b50% 的时间，但这种类型的结果在实践中基本上是无用的，因为你只会有N=2，并且统计推断不适用（即，每个组内没有可用于量化抽样不确定性的方差）。

您缺少的是比较两个感兴趣的汇总统计数据以确定它们的采样属性，这通常是统计推断的主题，并且至少需要最少数量的数据点来量化某种形式的采样不确定性。

就目前而言，这通常是独立于标准的t- 测试设置。因此，您可以比较许多事物，例如第一组平均值的频率高于第二组平均值的频率，t-test 给出一个p<α结果（或类似地，多久|t|比率大于截止值），依此类推。

例如，如果n=20对于每个组，然后在总体中的分布是相等的

runs <- 10000
n <- 20
alpha <- .05
result.p <- result.mean <- numeric(runs)

for(i in 1:runs){
  a <- rnorm(n, mean = 0) 
  b <- rnorm(n, mean = 0)
  result.mean[i] <- mean(a) > mean(b)
  result.p[i] <- t.test(a, b)$p.value
} 
mc.p.value <- mean(result.p < alpha)
mc.a_gt_b.value <- mean(result.mean)

使用其他参数，例如将第一组均值更改为 1，将改变模拟的性质（将其从 I 类错误模拟，就像现在一样，改为功率研究）。

其它你可能感兴趣的问题
上一篇贝叶斯与频率论：机器学习的实际差异 下一篇深度学习：同时使用 L2 和 Dropout 正则化？