我计划做一个实验来确定一个二元变量（值为1或0）的频率。
每天有10,000个新事件发生。
每天，我可以从新的一万人中随机抽出100人，看看他们的结果（不是1就是0
如何用这个数据估计人口中1和0的频率？
R中是否有一个软件包可以对这些数据进行离散概率分布拟合？
假设你有一个大小为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个。人口中总事件的最大似然估计值是简单的70100×10000=7000。这个估计的标准误差由以下公式给出
sqrt((N-n)/N * N^2 * var(sam$event)/n)
#[1] 473.71
95%的置信区间是[6101 - 7898]，涵盖了真实的人口总数6500人。但20天中有1天你很可能会得到一个坏的样本。
R包？对于这个实验来说，其实没有必要。对于更复杂的抽样设计，我只能想到的是 调查 包，但可能还有其他的。
现在，如果你重复这样做，比如说10天，你会得到每一天的估计值。根据频繁主义统计学家的说法，对总量的估计是总量×N n，对SE的估计也是以类似的方式计算。例如，假设你在连续十天从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
(0,1)事件的向量可以用以下方法构造： rep.
events <- rep(rep(c(0,1), each=10), times=c(events0, events1))
让我们将n和N分别定义为你的十天样本和十天人群中的事件数量。
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)相同。一个有经验的统计学家(频数学家或贝叶斯学家)会认识到贝塔分布为 共轭 的二项分布。因此，贝叶斯人将使用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% 而95%的概率是真实比例在7.84%和11.48%之间。以總人口推算，我們可以說，在10萬個事件中，9.58%或9,581個事件是正面的。正如我所说，这与频繁主义的方法非常相似。
元分析
现在，这两种方法实际上是把所有的十天都合并成一个大样本，估计阳性事件的比例，或者阳性事件的总数，在整个人群中的比例。更直观的做法可能是将每一天的结果以更合适的方式结合起来，基于权重，比如在荟萃分析中的做法。
如果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%的事件是阳性的，这与上述两种粗略的方法没有太大区别。由于十天样本的偏斜性，我们得到的估计值较小（也许更准确）。