R 中随机子集的频率 table

Question

假设我有一个 2006 年在美国出生的婴儿最流行的十个名字的频率图表：

myfreq <- c(24835, 22630, 22313, 21398, 20504, 20326, 20054, 19711, 19672, 19400)
names(myfreq) <- c("Jacob", "Michael", "Joshua", "Emily", "Ethan", "Matthew", "Daniel", "Andrew", "Christopher", "Anthony")

> myfreq
      Jacob     Michael      Joshua       Emily       Ethan     Matthew      Daniel 
      24835       22630       22313       21398       20504       20326       20054 
     Andrew Christopher     Anthony 
      19711       19672       19400

现在考虑 2006 年在美国出生的 210,843 个具有这些名字的婴儿。这个集合有 2^210843 个子集。我想要婴儿的 随机子集 的婴儿名字频率图表，每个子集的可能性均等。我的代码如下：

subfreq <- sapply(myfreq, function(k) sum(rbinom(k, 1, 0.5)))

这是在做我想让它做的事吗？有什么方法可以提高性能吗？它将处于具有数百万次迭代的循环中，并且 rbinom 函数似乎非常慢；我想知道对于 p=1/2 的二项式分布的这种特殊情况，R 中是否有更快的函数。感谢您的帮助。

Answer 1

不确定您的意思是不是想用引导程序模拟平局，但如果那是您想要的，我会尝试使用 data.table 的以下方法。单次开奖：

library(data.table)

# Example data:
dat.namefreqs <- data.table(name=LETTERS, count=sample(1e4, size=26))

# Format:
name count
   A  7466
   B 10000
   C  8897
   D  6833
   E  8614
   F  8128
   G  1837
   H  9349
   I  7798
   J  1158
   K  1707
   L  3368
   M  1019
   N   795
   O  1840
   P  4476
   Q  5345
   R   247
   S  5430
   T  9879
   U  1328
   V  4530
   W  6865
   X  6693
   Y  2186
   Z  1754

# Total all individuals
N.tot <- sum(dat.namefreqs$count)

# Repeat each name * its frequency
dat.expanded <- dat.namefreqs[rep(1:.N, count)]

# For a single random draw,
# Create a vector of binomial draws of 1s and 0s from rbinom, size = N.tot
# Use that as a true/false vector to extract names, and aggregate counts by name

dat.expanded[which(rbinom(N.tot, 1, 0.5)==1)][, .N, by=name]

单次抽取的示例输出：

    name    N
 1:    A 1339
 2:    B 1851
 3:    C 2898
 4:    D 4548
 5:    E 1066
 6:    F 4421
 7:    G 4754
 8:    H 3337
 9:    I 3144
10:    J  286
11:    K 1065
12:    L  880
13:    M 3435
14:    N 1942
15:    O 3851
16:    P 2471
17:    Q 3549
18:    R 4933
19:    S 1911
20:    T 3799
21:    U 4632
22:    V 1092
23:    W 3229
24:    X  631
25:    Y 1321
26:    Z 1883

并通过使用 foreach 进行引导来重复：我的机器在 17 秒内在单核上运行 ~1000 个引导程序，上面 table（136654 行，比你的一半多一点）

library(foreach)

dat.namefreqs <- data.table(name=LETTERS, count=sample(1e4, size=26))

N.tot <- sum(dat.namefreqs$count)

dat.expanded <- dat.namefreqs[rep(1:.N, count)]

results <- foreach(n=1:1000, .combine="rbind") %do% {
    dat <- dat.expanded[which(rbinom(N.tot, 1, 0.5)==1)][, .N, by=name]
    dat[, bootstrap := n]
    return(dat[])
}

> results
       name    N bootstrap
    1:    A 1339         1
    2:    B 1851         1
    3:    C 2898         1
    4:    D 4548         1
    5:    E 1066         1
   ---
25996:    V 1055      1000
25997:    W 3234      1000
25998:    X  636      1000
25999:    Y 1315      1000
26000:    Z 1895      1000

Answer 2

无法准确完成。您无法构造所有可能的子集，因此请忘记该方法。

如果你懂一些数学，可以大致完成。

首先你需要样本量为 n 的概率，这是（在 R 中）天真的：

choose(N, n)/2^N

对于中等 N 和 n（例如 N=1050 和 n=525），这将被打破。所以你可以尝试对数，经过一些工作你得到（其中 lgamma 是 gamma 函数的对数，n+1 处的 gamma 函数与 n 相同！）由下式给出的概率：

exp(lgamma(N+1) - lgamma(n+1) - lgamma(N-n+1) - N*log(2))

为了将所有概率放入一个向量中，我们可以将其封装到一个函数中：

pmf <- function(N,n) {
  exp(lgamma(N+1) - lgamma(n+1) - lgamma(N-n+1) - N*log(2))
}

N <- sum(myfreq)
probs <- sapply(0:N, function(n) pmf(N,n))

请注意，大多数样本量的概率为 0（大约）。现在，对于 select 您的样本，您将首先根据 probs 中的概率选择一个样本大小，然后从姓名总体中选择该大小的样本。我们需要首先根据您提供的频率制作该人口。

mypop <- rep(mynames, myfreq)

样本本身：

sample(mypop, sample(0:N, 1, prob = probs))

要复制很多次：

k <- 100
samps <- replicate(k, sample(mypop, sample(0:N, 1, prob = probs)))

samps 是随机 selected 大小的样本列表。

请注意，只有 non-zero 概率被 selected 的样本大小是：

range(which(probs > 0))
#> 96603 114242

所以你的样本的属性不会像你想象的那么有趣。它们将非常接近婴儿名字的人口分布。从一开始就制造婴儿肯定更有趣。

R 中随机子集的频率 table

Frequency table for a random subset in R

r

frequency-distribution