# 设置参数

posts <- 2     # 发帖数

n1 <- 2          # 设置生魂数量 = n1（计算总生魂数恰好为n1的概率）

n2 <- 0          # 设置生魂数量 ≤ n2（计算生魂数 ≤ n2的累积概率）

p1 <- 0.02       # 2%魂勋章的概率

p2 <- 0.01       # 1%魂勋章的概率

n1_max <- 12          # 2%魂勋章的数量

n2_max <- 7         # 1%魂勋章的数量

dist1 <- function(k) dbinom(k, size = n1_max, prob = p1)

dist2 <- function(k) dbinom(k, size = n2_max, prob = p2)

total_dist <- function(m) {

  sum(sapply(0:n1_max, function(k1) dist1(k1) * dist2(m - k1) * (m - k1 >= 0)))

}

# 计算总生魂数恰好为n1和生魂数≤n2的概率

calc_probabilities <- function(n1, n2) {

  # PGF of 单次发帖生魂数

  pgf <- sapply(0:(n1_max + n2_max), total_dist)

  

  # 卷积

  conv_pgf <- Reduce(function(x, y) convolve(x, rev(y), type = "o"), rep(list(pgf), posts))

  

  # 获取恰好为n1的概率

  exact_prob_n1 <- conv_pgf[n1 + 1]  # 因为R的索引从1开始，因此位置是n1+1

  

  # 获取生魂数≤n2的累积概率

  cumulative_prob_n2 <- sum(conv_pgf[1:(n2 + 1)])  # 累积概率从0到n2

  return(list(exact_prob_n1 = exact_prob_n1, cumulative_prob_n2 = cumulative_prob_n2))

}

probabilities <- calc_probabilities(n1, n2)

cat("发了", posts, "帖 总生魂数=", n1, "的概率为:", probabilities$exact_prob_n1, "\n")

cat("发了", posts, "帖 总生魂数≤", n2, "的累积概率为:", probabilities$cumulative_prob_n2, "\n")