背景
问题中的代码尝试在伯努利试验上运行蒙特卡罗模拟，以使用 Wald 置信区间计算覆盖率。代码中的一个问题是，许多计算是根据单个观察而不是成功和失败的总和来执行的。R 主要是一个矢量处理器，代码不会将单个观察结果汇总为成功和失败的计数，以计算 Wald 置信区间。
这会导致代码始终为原始帖子中测试的样本大小的 p 值高于 0.01 的覆盖率生成 0。我们使用原始帖子中的代码来隔离将错误引入算法的位置。
m我们设置一个种子，为、n和赋值p，并尝试生成 10,000 个大小为 的伯努利试验n。
set.seed(95014)
m<-10000
n<-5
p<-0.01
x <- rbinom(m,size=1,prob = p)
此时x是一个包含 10,000 个 true = 1、false = 0 值的向量。
> table(x)
x
0 1
9913 87
但是，x不是5 次伯努利试验的 10,000 次样本运行。鉴于这一事实，原始代码中算法的所有后续处理都将是不正确的。
下一行代码计算 的值p.hat。这应该是样本中 5 个元素聚合的单个值，而不是 10,000 个元素的向量，除非 x 中的每个元素都代表 5 个元素的样本。
p.hat <- x/n
table(p.hat)
> table(p.hat)
p.hat
0 0.2
9913 87
的准确计算p.hat，将向量视为一个样本如下：
> p.hat <- sum(x)/length(x)
> p.hat
[1] 0.0087
...这非常接近我们之前在代码中指定的总体 p 值 0.01，但仍不代表样本大小为 5 的 10,000 次试验。相反，p.hat如上定义的那样，代表样本大小为 10,000 的一次伯努利试验。
修复代码的两个小改动
在为伯努利试验独立开发了蒙特卡洛模拟器后（详见下文），很明显，通过一些调整，我们可以修复原始帖子中的代码，使其产生有效的结果。
首先，我们在第一个参数中m乘以，因此产生的试验次数是样本大小的 10,000 倍。我们还将结果转换为具有 10,000 行和列的矩阵。nrbinom()n
其次，我们使用rowSums()将试验与成功计数相加，并将 10,000 个元素的结果向量除以，在给定样本大小的情况下n为 生成正确的值。p.hat更正后，其余代码将p.hat按原计划工作。
f3 <- function(n,probs) {
res1 <- lapply(n, function(i) {
setNames(lapply(probs, function(p) {
m<-10000
n<-i
p<-p
# make number of trials m*n, and store
# as a matrix of 10,000 rows * n columns
x <- matrix(rbinom(m*n,size=1,prob = p),nrow=10000,ncol=i)
# p.hat is simply rowSums(x) divided by n
p.hat <- rowSums(x)/n
lower.Wald <- p.hat - 1.96 * sqrt(p.hat*(1-p.hat)/n)
upper.Wald <- p.hat + 1.96 * sqrt(p.hat*(1-p.hat)/n)
p.in.CI <- (lower.Wald <p) & ( p < upper.Wald )
covprob1<- mean(p.in.CI)
covprob1
}),paste0("p=",probs))
})
names(res1) <- paste0("n=",n)
res1
}
f3(n=c(10,15,20,25,30,50,100,150,200),probs = c(0.01,0.4, 0.8))
...和输出：
> f3(n=c(10,15,20,25,30,50,100,150,200),probs = c(0.01,0.4, 0.8))
$`n=10`
$`n=10`$`p=0.01`
[1] 0.0983
$`n=10`$`p=0.4`
[1] 0.9016
$`n=10`$`p=0.8`
[1] 0.8881
$`n=15`
$`n=15`$`p=0.01`
[1] 0.1387
$`n=15`$`p=0.4`
[1] 0.9325
$`n=15`$`p=0.8`
[1] 0.8137
$`n=20`
$`n=20`$`p=0.01`
[1] 0.1836
$`n=20`$`p=0.4`
[1] 0.9303
$`n=20`$`p=0.8`
[1] 0.9163
$`n=25`
$`n=25`$`p=0.01`
[1] 0.2276
$`n=25`$`p=0.4`
[1] 0.94
$`n=25`$`p=0.8`
[1] 0.8852
$`n=30`
$`n=30`$`p=0.01`
[1] 0.2644
$`n=30`$`p=0.4`
[1] 0.9335
$`n=30`$`p=0.8`
[1] 0.9474
$`n=50`
$`n=50`$`p=0.01`
[1] 0.3926
$`n=50`$`p=0.4`
[1] 0.9421
$`n=50`$`p=0.8`
[1] 0.9371
$`n=100`
$`n=100`$`p=0.01`
[1] 0.6313
$`n=100`$`p=0.4`
[1] 0.9495
$`n=100`$`p=0.8`
[1] 0.9311
这些结果看起来更像我们对模拟的期望：在低 p 值/小样本量下覆盖率较差，对于给定的 p 值，覆盖率随着样本量的增加而提高。
从头开始：一个 p 值/样本大小的基本模拟器
在这里，我们开发了一个解决方案，它以一组基本构建块迭代构建：一个 p 值、一个样本大小和 95% 置信区间。该模拟器还跟踪参数，因此我们可以将多个模拟的结果组合成易于阅读和解释的数据帧。
首先，我们创建了一个模拟器，用于测试从具有给定概率值的伯努利分布中抽取的 10,000 个大小样本。它聚合成功和失败，然后计算 Wald 置信区间，并生成输出数据框。出于模拟的目的，我们传递给模拟器的 p 值代表“真实”总体概率值。我们将看到模拟在置信区间中包含总体 p 值的频率。
我们设置参数来表示真实的总体 p 值为 0.5、样本量为 5、z 值为 1.96，表示 95% 的置信区间。我们为这些常量创建了函数参数，因此我们可以在后续代码中改变它们。我们还用于set.seed()使结果具有可重复性。
set.seed(90125)
simulationList <- lapply(1:10000,function(x,p_value,sample_size,z_val){
trial <- x
successes <- sum(rbinom(sample_size,size=1,prob = p_value))
observed_p <- successes / sample_size
z_value <- z_val
lower.Wald <- observed_p - z_value * sqrt(observed_p*(1-observed_p)/sample_size)
upper.Wald <- observed_p + z_value * sqrt(observed_p*(1-observed_p)/sample_size)
data.frame(trial,p_value,observed_p,z_value,lower.Wald,upper.Wald)
},0.5,5,1.96)
此代码与原始问题中的代码之间的一个关键区别在于，我们从 5 个样本中抽取样本rbinom()并立即将真实值的数量相加以计算成功的数量。这使我们可以计算observed_p为successes / sample_size。p.hat现在我们有了原始问题中所称内容的经验生成版本。
结果列表包括一个总结每个试验结果的数据框。
我们将数据帧列表组合成一个数据帧do.call()
simulation_df <- do.call(rbind,simulationList)
此时simulation_df是一个包含 10,000 行和 6 列的数据框。每行代表一次sample_size伯努利试验模拟的结果。我们将打印前几行来说明数据框的内容。
> dim(simulation_df)
[1] 10000 6
> head(simulation_df)
trial p_value observed_p z_value lower.Wald upper.Wald
1 1 0.5 0.6 1.96 0.17058551 1.0294145
2 2 0.5 0.2 1.96 -0.15061546 0.5506155
3 3 0.5 0.6 1.96 0.17058551 1.0294145
4 4 0.5 0.2 1.96 -0.15061546 0.5506155
5 5 0.5 0.2 1.96 -0.15061546 0.5506155
6 6 0.5 0.4 1.96 -0.02941449 0.8294145
>
请注意这些observed_p值是如何以 0.2 为增量的不同值。这是因为当样本大小为 5 时，每个样本中 TRUE 值的数量可以在 0 到 5 之间变化。直方图observed_p可以清楚地说明这一点。
即使样本大小为 5，我们也可以在直方图中看到二项分布的形状。
p_value接下来，我们通过对总体 p 值（表示为）在 Wald 置信区间内的行求和来计算覆盖率。
# calculate coverage: % of simulations where population p-value is
# within Wald confidence limits generated via simulation
sum(simulation_df$p_value > simulation_df$lower.Wald &
simulation_df$p_value < simulation_df$upper.Wald) / 10000 * 100
> sum(simulation_df$p_value > simulation_df$lower.Wald &
+ simulation_df$p_value < simulation_df$upper.Wald) / 10000 * 100
[1] 93.54
考虑到我们计算了 95% 的置信区间，93.54% 的覆盖率是合理的模拟。我们将此解释为 93.5% 的样本生成了 Wald 置信区间，其中包括 0.5 的总体 p 值。
因此，我们得出结论，我们的模拟器似乎正在生成有效的结果。我们将在此基本设计的基础上执行具有多个 p 值和样本大小的模拟。
模拟给定样本大小的多个 p 值
接下来，我们将改变概率值以查看 5 个观测值的 10,000 个样本的覆盖百分比。由于Sauro 和 Lewis, 2005等统计文献告诉我们，Wald 置信区间对于非常低和非常高的 p 值的覆盖率都很差，因此我们添加了一个参数来计算调整后的 Wald 分数。我们暂时将这个参数设置FALSE为。
p_val_simulations <- lapply(c(0.01,0.1,0.4,.5,.8),function(p_val){
aSim <- lapply(1:10000,function(x,p_value,sample_size,z_val,adjWald){
trial <- x
successes <- sum(rbinom(sample_size,size=1,prob = p_value))
if(adjWald){
successes <- successes + 2
sample_size <- sample_size + 4
}
observed_p <- sum(successes) / (sample_size)
z_value <- z_val
lower.Wald <- observed_p - z_value * sqrt(observed_p*(1-observed_p)/sample_size)
upper.Wald <- observed_p + z_value * sqrt(observed_p*(1-observed_p)/sample_size)
data.frame(trial,p_value,sample_size,observed_p,z_value,adjWald,lower.Wald,upper.Wald)
},p_val,5,1.96,FALSE)
# bind results to 1 data frame & return
do.call(rbind,aSim)
})
结果列表p_val_simulations包含一个数据框，用于模拟运行的每个 p 值。
我们组合这些数据框并计算覆盖率百分比如下。
do.call(rbind,lapply(p_val_simulations,function(x){
p_value <- min(x$p_value)
adjWald <- as.logical(min(x$adjWald))
sample_size <- min(x$sample_size) - (as.integer(adjWald) * 4)
coverage_pct <- (sum(x$p_value > x$lower.Wald &
x$p_value < x$upper.Wald) / 10000)*100
data.frame(p_value,sample_size,adjWald,coverage_pct)
}))
正如预期的那样，我们离 p 值 0.5 越远，覆盖率就越差。
p_value sample_size adjWald coverage_pct
1 0.01 5 FALSE 4.53
2 0.10 5 FALSE 40.23
3 0.40 5 FALSE 83.49
4 0.50 5 FALSE 94.19
5 0.80 5 FALSE 66.35
然而，当我们用 重新运行模拟时adjWald = TRUE，我们得到以下结果。
p_value sample_size adjWald coverage_pct
1 0.01 5 TRUE 95.47
2 0.10 5 TRUE 91.65
3 0.40 5 TRUE 98.95
4 0.50 5 TRUE 94.19
5 0.80 5 TRUE 94.31
这些要好得多，特别是对于接近分布末端的 p 值。
剩下的最后一项任务是修改代码，以便在不同级别的样本量下执行蒙特卡罗模拟。在继续之前，我们计算到目前为止我们开发的代码的运行时间。
system.time()告诉我们，在配备 2.5 Ghz Intel i-7 处理器的 MacBook Pro 15 上运行 10,000 次伯努利试验的 5 次不同蒙特卡罗模拟（样本大小为 5）的代码大约需要 38 秒。因此，我们预计下一次模拟将需要几分钟才能运行。
改变 p 值和样本量
我们添加了另一个级别lapply()来解释样本量的变化。我们还将adjWald参数设置为，FALSE以便我们可以看到基本 Wald 置信区间在 p = 0.01 和 0.10 时的表现。
set.seed(95014)
system.time(sample_simulations <- lapply(c(10, 15, 20, 25, 30, 50,100, 150, 200),function(s_size){
lapply(c(0.01,0.1,0.8),function(p_val){
aSim <- lapply(1:10000,function(x,p_value,sample_size,z_val,adjWald){
trial <- x
successes <- sum(rbinom(sample_size,size=1,prob = p_value))
if(adjWald){
successes <- successes + 2
sample_size <- sample_size + 4
}
observed_p <- sum(successes) / (sample_size)
z_value <- z_val
lower.Wald <- observed_p - z_value * sqrt(observed_p*(1-observed_p)/sample_size)
upper.Wald <- observed_p + z_value * sqrt(observed_p*(1-observed_p)/sample_size)
data.frame(trial,p_value,sample_size,observed_p,z_value,adjWald,lower.Wald,upper.Wald)
},p_val,s_size,1.96,FALSE)
# bind results to 1 data frame & return
do.call(rbind,aSim)
})
}))
MacBook Pro 上的经过时间为 217.47 秒，或约 3.6 分钟。鉴于我们运行了 27 次不同的蒙特卡洛模拟，代码每 8.05 秒完成一次模拟。
最后一步是处理列表列表以创建总结分析的输出数据框。我们聚合内容，将行组合成数据框，然后绑定数据框的结果列表。
summarizedSimulations <- lapply(sample_simulations,function(y){
do.call(rbind,lapply(y,function(x){
p_value <- min(x$p_value)
adjWald <- as.logical(min(x$adjWald))
sample_size <- min(x$sample_size) - (as.integer(adjWald) * 4)
coverage_pct <- (sum(x$p_value > x$lower.Wald &
x$p_value < x$upper.Wald) / 10000)*100
data.frame(p_value,sample_size,adjWald,coverage_pct)
}))
})
results <- do.call(rbind,summarizedSimulations)
最后一步，我们按 p 值对数据进行排序，以查看覆盖率如何随着样本量的增加而提高。
results[order(results$p_value,results$sample_size),]
...和输出：
> results[order(results$p_value,results$sample_size),]
p_value sample_size adjWald coverage_pct
1 0.01 10 FALSE 9.40
4 0.01 15 FALSE 14.31
7 0.01 20 FALSE 17.78
10 0.01 25 FALSE 21.40
13 0.01 30 FALSE 25.62
16 0.01 50 FALSE 39.65
19 0.01 100 FALSE 63.67
22 0.01 150 FALSE 77.94
25 0.01 200 FALSE 86.47
2 0.10 10 FALSE 64.25
5 0.10 15 FALSE 78.89
8 0.10 20 FALSE 87.26
11 0.10 25 FALSE 92.10
14 0.10 30 FALSE 81.34
17 0.10 50 FALSE 88.14
20 0.10 100 FALSE 93.28
23 0.10 150 FALSE 92.79
26 0.10 200 FALSE 92.69
3 0.80 10 FALSE 88.26
6 0.80 15 FALSE 81.33
9 0.80 20 FALSE 91.88
12 0.80 25 FALSE 88.38
15 0.80 30 FALSE 94.67
18 0.80 50 FALSE 93.44
21 0.80 100 FALSE 92.96
24 0.80 150 FALSE 94.48
27 0.80 200 FALSE 93.98
>
解释结果
蒙特卡洛模拟表明，即使样本量为 200，Wald 置信区间在 p 值为 0.01 时覆盖率也很差。覆盖率在 p 值为 0.10 时提高，在样本量为 25 及以上的模拟中，除了一个模拟之外的所有模拟超过 90%。对于 0.80 的 p 值，覆盖率甚至更好，其中除了一个样本大小超过 15 之外，所有样本量都超过了 91% 的覆盖率。
当我们计算调整后的 Wald 置信区间时，覆盖率会进一步提高，尤其是在 p 值较低的情况下。
results[order(results$p_value,results$sample_size),]
p_value sample_size adjWald coverage_pct
1 0.01 10 TRUE 99.75
4 0.01 15 TRUE 98.82
7 0.01 20 TRUE 98.30
10 0.01 25 TRUE 97.72
13 0.01 30 TRUE 99.71
16 0.01 50 TRUE 98.48
19 0.01 100 TRUE 98.25
22 0.01 150 TRUE 98.05
25 0.01 200 TRUE 98.34
2 0.10 10 TRUE 93.33
5 0.10 15 TRUE 94.53
8 0.10 20 TRUE 95.61
11 0.10 25 TRUE 96.72
14 0.10 30 TRUE 96.96
17 0.10 50 TRUE 97.28
20 0.10 100 TRUE 95.06
23 0.10 150 TRUE 96.15
26 0.10 200 TRUE 95.44
3 0.80 10 TRUE 97.06
6 0.80 15 TRUE 98.10
9 0.80 20 TRUE 95.57
12 0.80 25 TRUE 94.88
15 0.80 30 TRUE 96.31
18 0.80 50 TRUE 95.05
21 0.80 100 TRUE 95.37
24 0.80 150 TRUE 94.62
27 0.80 200 TRUE 95.96
调整后的 Wald 置信区间在 p 值和样本大小范围内始终提供更好的覆盖率，在 27 次模拟中平均覆盖率为 96.72%。这与表明调整后的 Wald 置信区间比未调整的 Wald 置信区间更保守的文献一致。
在这一点上，我们有一个有效的蒙特卡罗模拟器，它可以为多个 p 值和样本大小产生有效的结果。我们现在可以查看代码以寻找优化其性能的机会。
优化解决方案
遵循Make it work, make it right, make it fast的旧编程格言，以迭代的方式解决方案帮助我开发出产生有效结果的解决方案。
了解如何使它正确不仅使我能够看到问题中发布的代码中的缺陷，而且还使我能够设想解决方案。该解决方案使用rbinom()一次，参数为m * n，将结果转换为matrix()，然后用于计算 p 值，这让我了解如何通过消除每次模拟rowSums()的数千次调用来优化自己的解决方案。rbinom()
重构性能
我们创建了一个函数 ，binomialSimulation()它通过一次调用 来生成伯努利试验和 Wald 置信区间rbinom()，而不管单个模拟中的试验次数如何。我们还汇总了结果，因此每次模拟都会生成一个数据框，其中包含一行描述测试结果的行。
set.seed(90125)
binomialSimulation <- function(trial_size,p_value,sample_size,z_value){
trials <- matrix(rbinom(trial_size * sample_size,size=1,prob = p_value),
nrow = trial_size,ncol = sample_size)
observed_p <- rowSums(trials) / sample_size
lower.Wald <- observed_p - z_value * sqrt(observed_p*(1-observed_p)/sample_size)
upper.Wald <- observed_p + z_value * sqrt(observed_p*(1-observed_p)/sample_size)
coverage_pct <- sum(p_value > lower.Wald &
p_value < upper.Wald) / 10000 * 100
data.frame(sample_size,p_value,avg_observed_p=mean(observed_p),coverage_pct)
}
我们以 0.5 的总体 p 值、5 的样本量、10,000 次试验和 95% 的置信区间运行该函数，并使用 跟踪执行时间system.time()。优化后的函数比本文前面描述的原始实现快 99.8%，后者运行时间约为 6.09 秒。
system.time(binomialSimulation(10000,0.5,5,1.96))
> system.time(binomialSimulation(10000,0.5,5,1.96))
user system elapsed
0.015 0.000 0.015
我们将跳过中间步骤并展示迭代开发解决方案的优化版本。
system.time(results <- do.call(rbind,lapply(c(5,10,15,20,25,50,100,250),
function(aSample_size,p_values) {
do.call(rbind,lapply(p_values,function(a,b,c,d){
binomialSimulation(p_value = a,
trial_size = b,
sample_size = aSample_size,
z_value = d)
},10000,5,1.96))
},c(0.1,0.4,0.8))))
正如预期的那样，消除数千个不必要的调用rbinom()从根本上提高了解决方案的性能。
user system elapsed
0.777 0.053 0.830
鉴于我们之前的解决方案在 217 秒内运行，优化版本的性能确实令人印象深刻。现在我们有了一个解决方案，它不仅可以生成准确的伯努利试验蒙特卡罗模拟，而且速度也很快。