R - 如何加速递归和双重求和

Question

由于这本质上是一个关于如何在 R 中有效执行计算的问题，因此我将从方程开始，然后在代码之后为那些认为有用或有趣的人提供对该问题的解释。

我在 R 中编写了一个脚本，使用以下函数生成值:

正如您所看到的，该函数是递归的并且涉及双重求和。它对于 15 或更低的小数字效果很好，但在 n 和 t 值较高时，执行时间会变得非常长。我需要能够对从 1 到 30 的每个 n 和 t 对执行计算。有没有办法编写一个不需要几个月执行的脚本？

我当前的脚本是:

explProb <- function(n,t) {
    prob <- 0

    #################################
    # FIRST PART - SINGLE SUMMATION
    #################################
    i <- 0
    if(t<=n) {
        i <- c(t:n)
    }
    prob = sum(choose(n,i[i>0])*((1/3)^(i[i>0]))*((2/3)^(n-i[i>0])))

    #################################
    # SECOND PART - DOUBLE SUMMATION
    #################################
    if(t >= 2) {
        for(k in 1:(t-1)) {
            j <- c(0:(k-1))
            prob = prob + sum(choose(n,n-k)*((1/6)^(j))*((1/6)^(k-j))*((2/3)^(n-k))*explProb(k-j,t-k))
        }
    }

    return(prob)
}
MAX_DICE = 30
MAX_THRESHOLD = 30
probabilities = matrix(0,MAX_DICE,MAX_THRESHOLD)

for(dice in 1:MAX_DICE) {
    for(threshold in 1:MAX_THRESHOLD) {
        #print(sprintf("DICE = %d : THRESH = %d", dice, threshold))
        probabilities[dice,threshold] = explProb(dice,threshold)
    }
}

我正在尝试编写一个脚本来为桌面角色扮演游戏(具体来说，Shadowrun 第五版)中特定类型的骰子滚动生成一组概率。这种掷骰子的类型称为“爆炸掷骰子”。如果您不熟悉这些卷在游戏中的工作原理，请让我简要解释一下。

每当你尝试完成一项任务时，你都会通过掷一些六面骰子来进行测试。您的目标是在掷骰子时获得预定的“命中”次数。 “命中”被定义为六面骰子上的 5 或 6。因此，举例来说，如果您的骰子池中有 5 个骰子，并且您掷出:1、3、3、5、6，那么您将获得 2 次命中。

在某些情况下，您可以重新滚动所有已滚动的 6，以便尝试获得更多命中。这称为“爆炸”滚动。 6 算作命中，但可以重新滚动以“爆炸”成更多命中。为了清楚起见，我将举一个简单的例子......

如果您掷 10 个骰子，结果为 1, 2, 2, 4, 5, 5, 6, 6, 6, 6，那么您在第一次掷骰子时就得到了 6 次...然而，这 4 个骰子掷出的 6 可以再次重新掷出。如果您掷骰子并得到 3、5、6、6，那么您还有 3 次命中，总共 9 次命中。但你现在可以重新掷出另外两个 6...等等...你继续重新掷出 6，将 5 和 6 添加到你的总命中中，并继续下去，直到掷出没有 6 的掷骰。

上面列出的函数通过输入“骰子数量”和“命中次数”(此处称为“阈值”)来生成这些概率。

n = 掷骰子的数量
t = 要达到的“点击”阈值

Answer 1

使用转移矩阵计算

如果我们有 n=10 个骰子，则 prob=2 的事件发生 0 到 10 的概率/6 可以在 R 中有效地计算为

dbinom(0:10,10,2/6)

由于允许您继续滚动直到失败，因此任意数量的最终命中都是可能的(分布的支持是[0,Inf))，尽管概率呈几何递减。由于需要建立机器精度的截止值并且存在审查阈值，因此递归数值解决方案是可行的。

由于重新掷骰子的数量较少，因此预先计算所有转移概率是有意义的。

X<-outer(0:10,0:10,function(x,size) dbinom(x,size,2/6))

其中第 j 列的第 i 行给出了 (i-1) 成功(命中)的概率，其中(j-1) 试验(掷骰子)。例如，6 次试验恰好 1 成功的概率位于 X[2,7]。

现在，如果您从 10 骰子开始，我们可以将其表示为向量

d<-c(rep(0,10),1) 

表明，在概率为 1 的情况下，我们在其他地方都有 10 个骰子，概率为 0。

一次掷骰后，活骰子数量的概率为X %*% d。 两次掷骰后，概率为 X %*% X %*% d。我们可以通过迭代计算任意次数掷骰后的活骰子状态概率。

T<-Reduce(function(dn,n) X %*% dn,1:11,d,accumulate=TRUE)

其中 T[1] 给出了第一次掷骰之前活骰子的概率，T[11] 给出了 11< 之前活骰子的概率th(在 10th 之后)。

这足以计算预期值，但对于累积总和的分布，我们需要跟踪状态中的其他信息。以下函数在每一步 reshape 状态矩阵，使第 i 行和 j 列的概率为 (i-1) 当前累积总数为 j-1 的活骰子。

step<-function(m) {
  idx<-arrayInd(seq_along(m),dim(m))
  idx[,2]<-rowSums(idx)-1
  i<-idx[nrow(idx),]
  m2<-matrix(0,i[1],i[2])
  m2[idx]<-m
  return(m2)
}

为了恢复累积总数的概率，我们使用以下便利函数对反对角线求和

conv<-function(m) 
  tapply(c(m),c(row(m)+col(m)-2),FUN=sum)

继续快速滚动的概率会减小，因此我在 40 处截止，最多显示 20，四舍五入到 4 位

round(conv(Reduce(function(mn,n) X %*% step(mn), 1:40, X %*% d))[1:21],4)

#>      0      1      2      3      4      5      6      7      8      9 
#> 0.0173 0.0578 0.1060 0.1413 0.1531 0.1429 0.1191 0.0907 0.0643 0.0428 
#> 
#>     10     11     12     13     14     15     16     17     18     19   
#> 0.0271 0.0164 0.0096 0.0054 0.0030 0.0016 0.0008 0.0004 0.0002 0.0001  

<小时/>

模拟计算

这也可以使用简单的模拟在合理的时间内以合理的精度进行计算。

我们使用 sample(1:6,n,replace=TRUE) 模拟一掷 n 个 6 面骰子，计算重新掷骰子的数量，并迭代直到没有可用的，一路上计算“点击”。

sim<-function(n) {
  k<-0
  while(n>0) {
    roll<-sample(1:6,n,replace=TRUE)
    n<-sum(roll>=5)
    k<-k+n
   }
   return(k)
}

现在我们可以简单地复制大量试验并制成表格

prop.table(table(replicate(100000,sim(10))))

#>     0      1      2      3      4      5      6      7      8      9  
#> 0.0170 0.0588 0.1053 0.1431 0.1518 0.1433 0.1187 0.0909 0.0657 0.0421 
#>
#>     10     11     12     13     14     15     16     17     18     19
#> 0.0252 0.0161 0.0102 0.0056 0.0030 0.0015 0.0008 0.0004 0.0002 0.0001

即使使用 30 骰子，这也是相当可行的(即使重复 100,000 次也只需几秒钟)。