title: R - How to Speed Up Recursion and Double Summation
body: <p>Since this is essentially a question about how to efficiently perform a computation in R, I will start with the equation and then provide an explanation for the problem after the code for those who would find it useful or interesting. </p>
  
  <p>I have written a script in R to generate values using the following function:</p>
  
  <p><a href="https://i.stack.imgur.com/vLTPX.png" rel="noreferrer"><img src="https://i.stack.imgur.com/vLTPX.png" alt="Success probability when dice explode"></a></p>
  
  <p>The function, as you can see, is recursive and involves double summation. It works well for small numbers around 15 or lower, but the execution time gets prohibitively long at higher values of <code>n</code> and <code>t</code>. I need to be able to perform the calculation for every <code>n</code> and <code>t</code> pair from 1 to 30. Is there a way to write a script that won't take months to execute? </p>
  
  <p>My current script is:</p>
  
  <pre><code>explProb &lt;- function(n,t) {
      prob &lt;- 0
  
      #################################
      # FIRST PART - SINGLE SUMMATION
      #################################
      i &lt;- 0
      if(t&lt;=n) {
          i &lt;- c(t:n)
      }
      prob = sum(choose(n,i[i&gt;0])*((1/3)^(i[i&gt;0]))*((2/3)^(n-i[i&gt;0])))
  
      #################################
      # SECOND PART - DOUBLE SUMMATION
      #################################
      if(t &gt;= 2) {
          for(k in 1:(t-1)) {
              j &lt;- 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)
      }
  }
  </code></pre>
  
  <p>I am trying to write a script to generate a set of probabilities for a particular type of dice roll in a tabletop roleplaying game (Shadowrun 5th Edition, to be specific). The type of dice roll is called an "Exploding Dice Roll". In case you are not familiar with how these rolls work in this game, let me briefly explain. </p>
  
  <p>Whenever you try to accomplish a task you make a test by rolling a number of six-sided dice. Your goal is to get a predetermined number "hits" when rolling those dice. A "hit" is defined as a 5 or 6 on a six-sided die. So, for example, if you have a dice pool of 5 dice, and you roll: 1, 3, 3, 5, 6  then you have gotten 2 hits. </p>
  
  <p>In some cases you are allowed to re-roll all of the 6's that were rolled in order to try and get MORE hits.This is called an "exploding" roll. The 6's counts as hits, but can be re-rolled to "explode" into even more hits. For clarification I'll give a quick example... </p>
  
  <p>If you roll 10 dice and get a result of 1, 2, 2, 4, 5, 5, 6, 6, 6, 6 then you have gotten 6 hits on the first roll... However, the 4 dice that rolled 6's can be re-rolled again. If you roll those dice and get 3, 5, 6, 6 then you have 3 more hits for a total of 9 hits. But you can now re-roll the two more sixes you got... etc... You keep re-rolling the sixes, adding the 5's and 6's to your total hits, and keep going until you get a roll with no sixes. </p>
  
  <p>The function listed above generates these probabilities taking an input of "# of dice" and "number of hits" (called a "threshold" here).</p>
  
  <p><code>n = # of Dice being rolled
  t = Threshold number of "hits" to be reached</code></p>