置信区间的覆盖概率
Ham*_*aid 2 r function confidence-interval dataframe bernoulli-probability
从 Bernoulli(p) 中，我想计算各种样本大小 (n= 10, 15, 20, 25, 30, 50, 100, 150, 200) 的覆盖概率，以及 p = 0.01 时的每个样本大小， 0.4 和 0.8。
这是我的尝试，但除了 p=0.01 之外，到处都显示 0
f3 <- function(n,probs) {
res1 <- lapply(n, function(i) {
setNames(lapply(probs, function(p) {
m<-10000
n<-i
p<-p
x <- rbinom(m,size=1,p=p)
p.hat <- 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))
背景
问题中的代码尝试对伯努利试验运行蒙特卡罗模拟，以使用 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% 的样本生成了包含总体 p 值 0.5 的 Wald 置信区间。
因此，我们得出结论，我们的模拟器似乎正在生成有效的结果。我们将在此基本设计的基础上使用多个 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%。p 值为 0.80 时，覆盖率甚至更好，其中除了一个样本大小高于 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 值和样本大小生成有效结果。我们现在可以检查代码以寻找优化其性能的机会。
优化解决方案
遵循古老的编程格言“让它工作、让它正确、让它快速”，以迭代的方式制定解决方案帮助我开发了一个能够产生有效结果的解决方案。
了解如何使其正确不仅使我能够看到问题中发布的代码中的缺陷，而且还使我能够设想一个解决方案。该解决方案使用rbinom()一次带有 参数m * n，将结果转换为 a matrix()，然后用于rowSums()计算 p 值，让我了解如何通过消除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)
}
我们以总体 p 值为 0.5、样本大小为 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 秒，优化版本的性能确实令人印象深刻。现在我们有了一个解决方案，不仅可以生成准确的伯努利试验蒙特卡洛模拟，而且速度也很快。