估计 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 天内您可能会得到一个坏样品。

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

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

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

这意味着 "negative" 或未发生的事件是：

events0 <- 100 - events1

可以使用 rep.

如下构造 (0,1) 事件的向量

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

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

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

您的十天样本中 "positive" 事件的数量是：

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

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

贝叶斯会使用贝叶斯法则并使用统一 (0,1) 作为合理的先验分布，用于 "positive" 事件的比例。日期间。 Unif(0,1) 与 Beta(1,1) 相同。经验丰富的统计学家（频率派或贝叶斯派）会认识到 beta 分布与二项分布共轭。因此，贝叶斯将使用 Beta(1+X, 1+N-X) 分布计算十天期间 "positive" 事件的比例，其中 X 是 "positive" 事件的总数（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

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

估计 R 中的概率分布

Estimate Probability distribution in R

r

distribution

chi-squared

sampling