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 5th Edition,具体)中为特定类型的骰子滚动生成一组概率.骰子卷的类型称为"爆炸骰子卷".如果你不熟悉这些游戏在这个游戏中是如何工作的,那么让我简单解释一下.

每当你尝试完成一项任务时,你都可以通过滚动一些六面骰子进行测试.你的目标是在掷骰子时获得预定数量的"命中率"."击中"被定义为六面骰子上的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次命中.但你现在可以重新推动你得到的另外两个......等等......你继续重新推动六人,将5和6加到你的总命中,并继续前进,直到你得到一个没有六分的滚动.

上面列出的函数生成这些概率,输入"骰子数"和"命中数"(这里称为"阈值").

n = # of Dice being rolled t = Threshold number of "hits" to be reached

Answer 1

用转移矩阵计算

如果我们有n=10骰子,那么概率0到10与事件的发生prob=2/6可以在R上有效地计算为

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

由于你被允许继续滚动直到失败,所以可能有任何数量的最终命中(分布的支持[0,Inf)),尽管概率几何递减.递归数值解决方案是可行的,因为需要建立机器精度的截止值以及检查器的阈值的存在.

由于重新滚动的骰子数量较少,因此预先计算所有转移概率是有意义的.

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

i-th列的-th行j给出试验(i-1)成功(命中)的概率(j-1)(掷骰子).例如,试验1成功的概率6位于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]给出了活骰子的概率,并在11th 之前(在10th之后)给出了活骰子的概率.

这足以计算预期值,但对于累积总和的分布,我们需要跟踪该州的其他信息.以下函数在每一步重塑一个状态矩阵,以便i第-row和j-th列具有(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  

用模拟计算

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

我们模拟一组n6面骰子sample(1:6,n,replace=TRUE),计算重新滚动的数量,并迭代直到没有可用,沿途计数"命中".

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次重复).