如何在R中用sample()采样多项分布并计算特定概率？

一、用sample()完成多项分布采样 #

多项分布本质是n次独立有放回的类别采样，每个类别对应固定概率。你之前失败大概率是没开启有放回采样，或者没正确统计类别数量。以下是具体实现：

1. 单次采样 #

# 定义类别、采样规模和概率
categories <- c("X1", "X2", "X3")
sample_size <- 10
prob_vec <- c(0.3, 0.4, 0.3)

# 执行有放回采样（必须设置replace=TRUE）
sample_res <- sample(x = categories, size = sample_size, replace = TRUE, prob = prob_vec)

# 统计各类别出现次数，得到X1,X2,X3的采样结果
counts <- table(sample_res)
counts

2. 多次批量采样 #

如果需要重复采样（比如做模拟），用replicate()包裹采样逻辑：

# 模拟1000次采样
n_sim <- 1000
sim_counts <- replicate(n_sim, {
  res <- sample(x = 1:3, size = 10, replace = TRUE, prob = prob_vec)
  # 用factor强制保留所有类别，避免某次采样缺失类别导致table长度不足
  table(factor(res, levels = 1:3))
})

# 转成数据框方便查看结果
sim_df <- as.data.frame(t(sim_counts))
colnames(sim_df) <- categories
head(sim_df)

二、计算P(X1=3,X2=4,X3=3)的概率 #

分精确计算和模拟估算两种方式：

1. 精确计算（直接用多项分布公式） #

R内置的dmultinom()可以直接计算多项分布的概率密度：

exact_prob <- dmultinom(x = c(3, 4, 3), prob = prob_vec)
exact_prob

手动计算的话，公式是：
$$P(X_1=3,X_2=4,X_3=3) = \frac{10!}{3!4!3!} \times 0.3^3 \times 0.4^4 \times 0.3^3$$
用R代码实现手动计算：

manual_prob <- (factorial(10)/(factorial(3)*factorial(4)*factorial(3))) * (0.3^3) * (0.4^4) * (0.3^3)
manual_prob

2. 用sample()模拟估算概率 #

通过大量重复采样，统计出现(3,4,3)组合的频率来近似概率：

n_sim <- 100000  # 模拟次数越多，结果越精确
target_count <- 0

for (i in 1:n_sim) {
  res <- sample(x = 1:3, size = 10, replace = TRUE, prob = prob_vec)
  cnt <- table(factor(res, levels = 1:3))
  if (cnt[1]==3 && cnt[2]==4 && cnt[3]==3) {
    target_count <- target_count + 1
  }
}

estimated_prob <- target_count / n_sim
estimated_prob

内容的提问来源于stack exchange，提问作者Shiyu Wang