估计 R 中的概率分布

Question

我正在计划一个实验来确定一个二元变量（值为 1 或 0）的频率。

每天都有 10,000 个新事件发生

每天，我都会从新的 10,000 个中随机抽取 100 个并查看它们的结果（1 或 0）

如何使用这些数据估计总体中 1 和 0 的频率？

R 中是否有可以将离散概率分布拟合到此数据的包？

Answer 1

假设您的人口规模为 N=10,000，其中一天发生了 6,500 个事件。

pop <- rep(c(0,1), times=c(3500, 6500))
table(pop)
#pop
#   0    1 
#3500 6500

现在假设您可以对这些 (0,1) 事件中的 100 个进行采样无需替换。

set.seed(123)
N <- 10000
n <- 100
sam <- data.frame(id=1:n, event=sample(pop, size=n), prob=n/N)

table(sam$event)
# 0  1 
#30 70

所以我们在 100 个中得到了 70 个。总体中事件总数的最大似然估计是 70/100 x 10,000 = 7,000。此估计的标准误差由

给出

sqrt((N-n)/N * N^2 * var(sam$event)/n)
#[1] 473.71

95% 的置信区间是 [6101 - 7898]，它涵盖了 6,500 的真实人口总数。但是 20 天中有 1 天可能会得到一个坏样本。

R 包？这个实验真的没有必要。对于更复杂的抽样设计，我只能想到survey包，但可能还有其他的。

---

现在，如果您反复执行此操作，例如 10 天，您将获得每天的估算值。根据一位常客统计学家的说法，对总数的估计将是总数 x N / n 和以类似方式计算的 SE 的估计。例如，假设您连续 10 天从 100 个样本中发现 3、4、5、11、6、8、14、8、17 和 19 个“阳性”事件：

events1 <- c(3, 4, 5, 11, 6, 8, 14, 8, 17, 19)

这意味着“负面”或不发生的事件是：

events0 <- 100 - events1

可以使用rep 构造如下的 (0,1) 事件向量。

events <- rep(rep(c(0,1), each=10), times=c(events0, events1))

让我们将 n 和 N 分别定义为您的 10 天样本和 10 天总体中的事件数。

n <- 100 * 10
N <- 10000 * 10

十天样本中“积极”事件的数量为：

sum(events1)
#[1] 95

十天人口中的MLE为：

(T <- sum(events1) * N / n)
[1] 9500

这个十天估计的标准误是：

SE <- sqrt((N-n)/N * N^2 * var(events)/n); SE
[1] 923.0409

95% CI：

T + c(-1,1) * 1.96*SE
[1]  7690.84 11309.16

贝叶斯可能会说每天都应该根据前一天的估计重新估计或更新，但我认为结果会非常相似。

---

贝叶斯将使用贝叶斯规则并使用 Uniform(0,1) 作为合理的先验分布，用于十天期间的“积极”事件的比例。 Unif(0,1) 与 Beta(1,1) 相同。有经验的统计学家（Frequentist 或 Bayesian）会认识到 beta 分布与二项分布共轭。因此，贝叶斯将使用 Beta(1+X, 1+N-X) 分布来表示十天期间“积极”事件的比例，其中 X 是“积极”事件的总数 (95)，N 是样本量（1000）。注意 Beta(alpha, beta) = alpha/(alpha+beta) 的平均值。

在 R 中：

n <- rep(100, 10)
events1 <- c(3, 4, 5, 11, 6, 8, 14, 8, 17, 19)

X <- sum(events1)
N <- sum(n)

pmean = (1+X)/(2+N); pmean
#[1] 0.09580838

CI = qbeta(c(.025,.975), 1+X, 1+N-X); CI # 95% credible interval
#[1] 0.07837295 0.11477134

因此，在十天的时间里，阳性事件的比例将是所有事件的 9.58%，真实比例在 7.84% 和 11.48% 之间的概率为 95%。外推到总人口，我们可以说 100,000 个事件或 9,581 个事件中有 9.58% 是积极的。正如我所说，这与频率论方法非常相似。

元分析

现在，这两种方法有效地将所有十天组合成一个大样本，并估计阳性事件在整个人群中的比例，或阳性事件的总数。基于权重以更合适的方式组合每天的结果可能更直观，例如在元分析中所做的。

如果 p[k] 是第 k 天的估计比例，se[k] 是它的标准误，那么组合估计由 p_hat = sum(p[k] * w[k]) / sum( w[k])，其中 w[k] = (1 / se[k])^2，标准误差为 1 / sqrt(sum(w[k])。

在 R 中：

N <- rep(100000, 10) # Population and 10 day period
n <- rep(100, 10) 

events1 <- c(3, 4, 5, 11, 6, 8, 14, 8, 17, 19)
events0 <- n - events1

p <- NULL; SE <- NULL; w <- NULL

for(k in seq_along(events1)){
  events <- c(rep(0, events0[k]), rep(1, events1[k]))
  p[k] <- sum(events1[k]) / n[k]
  SE[k] <- sqrt((N[k]-n[k]) / N[k] * var(events)/n[k])
  w[k] <- 1 / (SE[k])^2
}

p_hat <- sum(p*w)/sum(w); p_hat
#[1] 0.06997464

SE_p <- 1 / sqrt(sum(w)); SE_p
#[1] 0.007943816

(p_hat + c(-1,1) * 1.96 * SE_p)
#[1] 0.05440476 0.08554452

因此，大约 7% 的事件将为阳性，置信区间为 95% (5.44% - 8.55%)，这与上面使用的两种粗略方法没有太大区别。由于十天样本的偏斜性质，我们得到了一个更小（也许更准确）的估计值。