Say we role $n$ identical, fair dice, each with $d$ sides (every side comes up with the same probability $\frac{1}{d}$). On each die, the sides are numbered from $1$ to $d$ with no repeating number, as you would expect. So an ordinary $d$ sided die pool.

Every dice in the outcome that shows a number equal or higher than the threshold number $t$ is said to show a hit. Every die that shows the maximum result of $d$ is rolled again, which we call "exploding". If the re-rolled dice show hits, the number of hits is added to the hit count. Dice that show the maximum after re-rolling are rolled again and their hits counted until none show a maximum result. Given the values of

$$ d\ ...\ \text{Number of sides on each die}\ \ d>0 $$ $$ n\ ...\ \text{Number of dies rolled}\ \ n\ge 0$$ $$ h\ ...\ \text{Number of hits, we want the probability for}$$ $$ t\ ...\ \text{Threshold value for a die to roll a hit}\ \ 0 < t \le d$$

what is the probability to get exactly exactly $h$ hits? Lets call it: $$p^\text{exploding}(d,n,t,h) = p_{d,n,t,h}$$ Can you derive a formula for this probability?

Example roll:

We roll 7 six-sided dice and count those as hits that show a 5 or a 6. In this example, $d=6$, $n=7$, $t=5$. The outcome of such a roll may be 6,5,1,2,3,6,1. That's three hits so far, but we have to roll the two sixes again (they explode). This time it's 6, 2. One more hit, and one more die to roll. We are at four hits at this point. The last die to be re-rolled shows 6 again, we re-roll it yet another time. On the last re-roll it shows a 4 - no more hits. That gives five hits in total and the roll is complete. So, for this roll $h=5$.

Simple case for just one die $n=1$:

If we roll only one die with the same threshold as above, so ($d=6$, $n=1$, $t=5$), the probabilities can be easily calculated:

$$ p_{6,1,5,0} = \frac{4}{6} \quad \text{(Probability for exactly 0 hits - roll 1-4 on the first roll, no explosion here)} $$ $$ p_{6,1,5,1} = \frac{1}{6} + \frac{1}{6} \cdot \frac{4}{6} \quad \text{(Probability for exactly 1 hit - roll either a 5 or a result of 1-4 after a 6)} $$ $$ p_{6,1,5,2} = \frac{1}{6} \cdot \frac{1}{6} + \frac{1}{6} \cdot \frac{1}{6} \cdot \frac{4}{6} \quad \text{(Probability for exactly 2 hits - either a 6 and 5 or two sixes and 1-4)} $$ $$ p_{d,1,t,h\ge 1} = \left(\frac{1}{d}\right)^{h-1}\frac{d-t}{d} + \left( \frac{1}{d} \right)^h \cdot \frac{t-1}{d} \quad \text{(Probability for exactly $h\ge 1$ hits - either $h-1$ maximum rolls and non-maximal success or $h$ maximum rolls and a non-success )} $$

Without Explosion:

For none-exploding dice the probability would just be binomially distributed:

$$ p^\text{non-exploding}_{d,n,t,h} = \binom{n}{h} \left( \frac{d-t+1}{d} \right)^h \left( 1 - \frac{d-t+1}{d} \right)^{n-h} $$

$$ E^\text{non-exploding}_{d,n,t} = n \frac{d-t+1}{d}; \qquad V^\text{non-exploding}_{d,n,t} = n \frac{(d-1)(d-t+1))}{d^2} $$

Where $E_{d,n,t}$ is the expected number of hits, and $V_{d,n,t}$ its variance.

Edit1: In the mean time I found Probability of rolling $n$ successes on an open-ended/exploding dice roll. However I'm afraid, I don't fully get the answer there. E.g. the author says $s = n^k + r$, which does not hold for his examples. Also I'm not sure how to get $s$, $k$ and $r$ from my input values stated above (which are $d$, $n$, $h$ and $s$).

Edit2: If one had the probability for $b$ successes via explosions, given that the initial role had $l$ successes prior to the explosions, one could just subtract all those probabilities for all values of $b$ from the value for the pure binomial distributions with $l$ successes and add the respective value to the pure binomial probability of $b+l$ successes. Just an idea. I suppose this should be something like a combination of geometric and binomial distribution.

Edit3: I accepted Brian Thug's excellent answer, giving the formula: $$ p^\text{exploding}_{d,n,t,h} = \frac{(t-1)^n}{d^{n+h}} \sum_{k=0}^{\max\{h, n\}} \binom{n}{k} \binom{n+h-k-1}{h-k} \left[ \frac{d(d-t)}{t-1} \right]^k $$

$$ E^\text{exploding}_{d,n,t} = n\frac{d+1-t}{d-1}; \qquad V^\text{exploding}_{d,n,t} = E_{d,n,t} - n\frac{(d-t)^2-1}{(d-1)^2} $$

Here is a graph from a simulation (html) that illustrates the whole thing:

Comparison between exploding and non-exploding dice pools

  • 2
    $\begingroup$ Can you explain why $p_{6,1,5,1} = \frac{2}{6} \cdot \frac{4}{6}$? It seems to me that there are two routes to get exactly $1$ hit: Route 1: Roll a $5$ (Probability of this route is $\frac{1}{6}$.) Route 2: Roll a $6$, which yields a hit but requires a reroll, and then roll something less than $5$ so there is no additional hit. (Probability of second route is $\frac{1}{6}\frac{4}{6}$.) The routes are exclusive, so the total probability of exactly $1$ hit is $\frac{1}{6}+\frac{1}{6}\frac{4}{6}=\frac{5}{18}$, but you say it’s $\frac{2}{9}$. I may be misunderstanding the problem statement. $\endgroup$ – Steve Kass Feb 8 '16 at 1:26
  • 1
    $\begingroup$ Now that I look further, I think $p_{6,1,5,2}$ is also wrong. There are two routes to getting two hits: Route 1: Roll a $6$, then a $5$. Route 2: Roll $6$, a second $6$, and then something between $1$ and $4$. Your value of $\frac{1}{36}$ only accounts for the first route. $\endgroup$ – Steve Kass Feb 8 '16 at 20:22
  • $\begingroup$ And $p_{6,1,5,h\ge2}$ was wrong as well. Either $h-1$ sixes and a five, the probability for which is $\frac{1}{6^h}$, or $h$ sixes and any non hit roll $\rightarrow \frac{1}{6^h} \frac{4}{6}$. $\endgroup$ – con-f-use Feb 8 '16 at 23:41
  • $\begingroup$ So that I understand it in terms that I'm familiar with (and maybe you are too): Each die is independent of the others. There is a "to-hit" score, namely $t$. If the die rolls $t$ or higher, it is counted as a hit; furthermore, a roll of $d$ (the maximum for the die) is counted as a "critical hit", whose bonus is to roll the die again. Is that all basically right? $\endgroup$ – Brian Tung Feb 10 '16 at 19:17
  • $\begingroup$ @con-f-use: By this point, reading the OP and all the other answers, I assumed that I had more or less understood the question correctly, but it's always nice to get confirmation. I believe my answer is essentially correct, but it'd be nice if someone else looked at it critically. $\endgroup$ – Brian Tung Feb 10 '16 at 23:59

ETA: OK, I think I've fixed the problem. Off-by-one error...

I think this can be done with generating functions. The generating function for a single die is given by

$$ F(z) = \frac{t-1}{d} + \frac{(d-t)z}{d} + \frac{zF(z)}{d} $$

We can interpret this as follows: The probability that there are no hits on the one die is $\frac{t-1}{d}$, so $F(z)$ has that as the coefficient for $z^0 = 1$. The probability that there is one hit and the die doesn't "explode" (repeat) is $\frac{d-t}{d}$, so $F(z)$ has that as the coefficient for $z^1 = z$. In the remaining $\frac{1}{d}$ of the cases, the die explodes and the situation is exactly as it was at the start, except that there is one hit already to our credit, which is why we have $zF(z)$: the $F(z)$ takes us back to the beginning, so to speak, and the multiplication by $z$ takes care of the existing hit.

This expression can be solved for $F(z)$ via simple algebra to yield

$$ F(z) = \frac{t-1+(d-t)z}{d-z} $$

whose $z^h$ coefficient gives the probability for $h$ hits. For example, for the simple case $n = 1, d = 20, t = 11$:

\begin{align} F(z) & = \frac{10+9z}{20-z} \\ & = \frac{10+9z}{20} \left(1+\frac{z}{20}+\frac{z^2}{20^2}+\cdots\right) \\ & = \left( \frac{1}{2} + \frac{9}{20}z \right) \left(1+\frac{z}{20}+\frac{z^2}{20^2}+\cdots\right) \\ \end{align}

and then we obtain the probability that there are $h$ hits from the $z^h$ coefficient of $F(z)$ as

$$ P(H = h) = \frac{1}{2\cdot20^h}+\frac{9}{20^h} = \frac{19}{2\cdot20^h} \qquad h > 0 $$

with the special case

$$ P(H = 0) = \frac{1}{2} $$

In general, we can obtain the expectation of the number of hits $\overline{H}$ as

$$ \overline{H} = F'(1) = \frac{d(d-t)+t-1}{(d-1)^2} = \frac{d+1-t}{d-1} $$

Now, for $n$ dice, we have

$$ [F(z)]^n = \left[ \frac{t-1+(d-t)z}{d-z} \right]^n $$

We can write this as $N(z)M(z)$, where

\begin{align} N(z) & = [t-1+(d-t)z]^n \\ & = \sum_{k=0}^n \binom{n}{k} (t-1)^{n-k}(d-t)^kz^k \end{align}


\begin{align} M(z) & = \left(\frac{1}{d-z}\right)^n \\ & = \frac{1}{d^n} \left( 1+\frac{z}{d}+\frac{z^2}{d^2}+\cdots \right)^n \\ & = \sum_{j=0}^\infty \binom{n+j-1}{j} \frac{z^j}{d^{n+j}} \end{align}

so we can obtain a closed form for $P(H = h)$ from the $z^h$ coefficient of $[F(z)]^n = N(z)M(z)$ as

\begin{align} P(H = h) & = \sum_{k=0}^{\max\{h, n\}} \binom{n}{k} \binom{n+h-k-1}{h-k} \frac{(t-1)^{n-k}(d-t)^k}{d^{n+h-k}} \\ & = \frac{(t-1)^n}{d^{n+h}} \sum_{k=0}^{\max\{h, n\}} \binom{n}{k} \binom{n+h-k-1}{h-k} \left[ \frac{d(d-t)}{t-1} \right]^k \end{align}

For example, for $n = 1, d = 6, t = 5$ (the example in the OP), the above expression yields

$$ P(H = h) = \frac{5}{3 \cdot 6^h} \qquad h > 0 $$

with the special case

$$ P(H = 0) = \frac{2}{3} $$

which coincides with the conclusions drawn in the comments to the OP.

The expectation for the number of hits could be obtained by evaluating $\frac{d}{dz} [F(z)]^n$ at $z = 1$, but owing to the linearity of expectation, it is obtained more straightforwardly as $n$ times the expected number of hits for one die, namely

$$ \overline{H} = \frac{n(d+1-t)}{d-1} $$

I think this all checks out, but some independent verification (or disproof, as appropriate) would be nice.

| cite | improve this answer | |
  • 1
    $\begingroup$ And now my colleague has come in and sown the seed of conceiving this question as a network congestion problem. $\endgroup$ – Brian Tung Feb 10 '16 at 23:10
  • $\begingroup$ I tested your results with a little python simulation. Will edit your post to include the result. I'd say, you're pretty spot on. Thank you very much. I think I learned a lot just going through your post. $\endgroup$ – con-f-use Feb 12 '16 at 14:09
  • $\begingroup$ Or not, as it seems that my edit was rejected for what ever reason. $\endgroup$ – con-f-use Feb 12 '16 at 14:47
  • $\begingroup$ Nice approach! (+1) $\endgroup$ – Markus Scheuer Feb 12 '16 at 15:14
  • $\begingroup$ @con-f-use: I think it would be OK to add your analysis as a separate answer. $\endgroup$ – Brian Tung Feb 12 '16 at 16:26

Each die produces a certain number of hits independently from the other dies. Denote by $p_r$ $(r\geq0)$ the probability that a given die produces exactly $r$ hits. Then $$p_0={t-1\over d},\qquad p_1={d-t \over d}+{1\over d} p_0\ ,$$ since the bonus roll (after a $d$) should not result in a hit. For $r>1$ hits we need $r-1$ rolls of a $d$ and exactly one more hit. It follows that we have $$p_r={1\over d^{r-1}}p_1={1\over d^r}(d-t+p_0)\qquad(r\geq1)\ .$$ Consider the generating function $$f(x):=\sum_{r=0}^\infty p_r x^r=p_0+(d-t+p_0){x/d\over 1-x/d}=p_0{1+\alpha x \over 1-\beta x}$$ with $$\alpha:={d-t\over d p_0},\qquad \beta:={1\over d}\ .$$ The probability $p_r^{(n)}$ of obtaining a total of $r$ hits with $n\geq1$ dice is the coefficient of $x^r$ in the expansion of $\bigl(f(x)\bigr)^n$. Now $$\bigl(f(x)\bigr)^n=p_0^n \sum_{j=0}^n {n\choose j}(\alpha x)^j\sum_{k=0}^\infty {n+k-1\choose k}(\beta x)^k\ .$$ When you collect terms you obtain $$p_r^{(n)}=p_0^n\sum_{j=0}^r C_j\>\alpha^j\beta^{r-j}$$ with certain combinatorial coefficients $C_j$. I don't know whether a simplification is possible.

The expected number $E^{(n)}$ of hits can be found more simply as follows: Let $E$ be the expected number $E$ of hits that you can land with a single die. Then $$E={t-1\over d}\cdot 0+{d-t\over d}\cdot 1+{1\over d}(1+E)\ .$$ The reasoning behind this is the following: If the first roll results is one of $\{1,\ldots,t-1\}$ there is no hit, and the die is dead. If the first roll results in one of $\{t,\ldots, d-1\}$ there is one hit, and the die is dead. If the first roll results in $d$ there is one hit, and in addition we are given the opportunity to realize $E$ more hits with this die.

It follows that $$E={d-t+1\over d-1}\ ,$$ from which we obtain $$E^{(n)}=n \>{d-t+1\over d-1}$$ when there are $n$ dice in the game.

| cite | improve this answer | |
  • 1
    $\begingroup$ +1. Essentially the same approach I took, it just took me a while to recognize they were the same. $\endgroup$ – Brian Tung Feb 12 '16 at 16:25

Fix $d$ and $t$, and abbreviate $p_{d,n,t,h}$ by $p_{n,h}$. Then:

$$p_{n,h}=\sum_{i=0}^{h}\sum_{j=0}^{i}q_{n,i,j}p_{j,h-i}$$ where $p_{0,0}=1$ and $q_{n,i,j}$ denotes the probability of $i$ hits and $j$ maximals by throwing $n$ dice.

For $q_{n,i,j}$ we find the expression:

$$q_{n,i,j}=\frac{n!}{j!\left(i-j\right)!\left(n-i\right)!}\times d^{-n}\left(h-t\right)^{i-j}\left(t-1\right)^{n-i}$$

| cite | improve this answer | |
  • $\begingroup$ So to get the result for one specific value of n, one first has to calculate it for every smaller value. Ouch! Didn't expect that. Also when $n>0$ there is a probability $>0$ for every value of $h$ even very large ones. That never stops. $\endgroup$ – con-f-use Feb 10 '16 at 19:02
  • $\begingroup$ But each of the maximals turns into some number of hits in the end, so this is only one step toward an eventual solution. $\endgroup$ – Brian Tung Feb 10 '16 at 21:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.