【派遣s1】用R语言计算先锋存活的概率

咸鱼鱼 · 发表于 4 天前

本帖最后由咸鱼鱼于 2025-4-25 13:24 编辑

先知：远征期间将有先知每日掷出生死骰子预测先锋的生死，累计结果。
         若掷出数量生(1)＞死(-1)，则先锋必将生还凯旋。
         若掷出数量生(1)＜死(-1)，则先锋必将英勇牺牲。
         若掷出数量生(1)= 死(-1)，则次日再掷1次定生死。
起止：先知预测有效期: 4月25日~5月7日

规则如上，可以简单翻译为
一个先知每天占卜先锋的生死，有三种情况，生、死、不明，概率相等。持续X日，可以暂时认为7日或者14日（毕竟7天后才选先锋，前面7天一直出死是不是大家就未战先怯了？）

当生的次数大于死的次数，则先锋生存。当死的次数大于生时，则先锋死亡。如果次数相同，则额外占卜若干次，直到生死次数不等。求先锋生存的概率

如果占卜仅有概率相等的两种情况（生、死），其他规则不变，求先锋生存的概率

此外在最开始的设计中，是d-1会出现-1、0、1三种情况的，我们也不知道最终结果是会保留三种情况，还是只保留生死，因此都必须计算

三种情况的代码

compute_prob_case1 <- function(X) {
# 生成所有可能的s和d组合
dat <- expand.grid(s = 0:X, d = 0:X)
dat <- dat[dat$s + dat$d <= X, ]
# 计算每个组合的概率
dat$prob <- apply(dat, 1, function(row) {
s <- row[['s']]
d <- row[['d']]
counts <- c(s, d, X - s - d)
dmultinom(counts, prob = rep(1/3, 3))
})
P_gt <- sum(dat$prob[dat$s > dat$d])
P_eq <- sum(dat$prob[dat$s == dat$d])
return(P_gt + 0.5 * P_eq)
}
X <- 7
result_case1 <- compute_prob_case1(X)
print(result_case1) # 输出：0.5

复制代码

两种情况的代码

compute_prob_case2 <- function(X) {
# 由于对称性，生存概率始终为0.5
return(0.5)
}
X <- 7
result_case2 <- compute_prob_case2(X)
print(result_case2) # 输出：0.5

复制代码

最终的结论
问题1分析：在每次占卜中，生、死、不明的概率均为1/3。当生死次数相等时，进入额外占卜阶段，此时每次占卜的生/死概率相等，故生存概率为0.5。由于生与死的初始概率对称，总体生存概率为0.5。

问题2分析：
每次占卜生或死的概率各为1/2。当生死次数相等时，额外占卜的生存概率为0.5。同样由于对称性，无论占卜天数X是奇数还是偶数，最终生存概率均为0.5。

两个问题的生存概率均为0.5，这可以通过对称性和概率的公平性解释，R代码验证了这一点。

匪夷所思的反直觉事件
为什么无论怎么样都是0.5呢，我觉得可以简单理解成抛硬币，你抛多少枚硬币，是不是无论每个面都是0.5的概率，那么在1次定胜负的概率是0.5的情况下，两次的情况是不是，正正，正反，反正，正正。还是0.5？。当你抛第三次的时候，在第二次的基础上
正正和反反，第三次硬币都不会改变第二次的结果。正反，反正则会基于第三次的结果发生改变，依旧是0.5

那么三种情况下，你就可以理解成不明的时候，是硬币立起来了，无需计算其中。还是只要看正反情况
所以，无论持续占卜多少天，先锋生存的概率都是0.5

死去的高数还在追我……

凯诺斯 · 发表于 4 天前

生存的概率都是0.5那就只能看运气了惹

毛茸茸兽兽 · 发表于 4 天前

一半的存活概率咩(´×ω×`)感觉要遭，而且还没那个复活的勋章的……说不定压根不会被抽到欸～那没事咯，喝茶ing

诺鸦cc · 发表于 4 天前

一半的概率）好残酷的LV7级用户随从远征.jpg大佬们加油啊）

野生阿努厨 · 发表于 4 天前

如果这个概率模型复杂一点就可以计算马尔可夫链了，不过看起来版主不打算折磨大家的大脑

好像还可以通过转经筒复活来着，要是一死就能复活那应该能提升不少概率

Apocatastasis · 发表于 4 天前

毕竟最后生死是统计生死数量，而不是全部都要roll生才能活。只要这个概率不涉及条件和状态的改变，那在这个赛博空间里，这个概率是不会随着次数递增受到影响的。【真实世界扔骰子就不一样了【当然，赛博伪随机数什么的就不讨论了
话说比起apply，这种不复杂的function我更喜欢直接向量计算。感觉这样简单一点
compute_prob_case1 <- function(X) {
  s_vec <- rep(0:X, each = X + 1)
  d_vec <- rep(0:X, times = X + 1)
  valid_idx <- s_vec + d_vec <= X

  s <- s_vec[valid_idx]
  d <- d_vec[valid_idx]
  r <- X - s - d

  probs <- dmultinom(x = cbind(s, d, r), prob = rep(1/3, 3))

  P_gt <- sum(probs[s > d])
  P_eq <- sum(probs[s == d])

  return(P_gt + 0.5 * P_eq)
}