# Roll and Keep in Anydice

I want to check out the probabilities of the Roll and Keep system (as known from L5R 1-4th Edition) in Anydice. We luckily have boundaries:

• We roll Skill + Stat + Modifier initial dice.
• 10 explodes; A roll of 10-10-1 is counted as a single die of 21!
• We keep (the best/worst) Ring dice of the roll.
• there are at best 10 dice rolled, any 2 dice more that could be rolled add 1 die kept.
• you can't ever keep more dice than you rolled.

For a roll of 10k3, this would be this code, but it fails to compute:

output [highest 3 of 10d [explode d10]] named "10k3"


How best to model this in a way that allows testing 1k1 to 10k5?

• A starting point for script development: Our Burning Wheel AnyDice question has some interesting information on how do to chained exploding dice correctly (it's hard): How can I model Burning Wheel probabilities in AnyDice? – Alex P Mar 20 '20 at 16:58
• Thank you for clarifying. – Wesley Obenshain Sep 1 '20 at 13:32

# Cut your losses accuracy

When working with many kinds of modeling there is a tradeoff between accuracy and speed. Anydice cuts us off on our speed, so we need to lose some accuracy. Anydice already truncates to exploding twice anyway and seeing as exploding twice on a d10 will only happen in 1 in a hundred throws the error should be fairly small.

For simplicity — rather than implementing a custom explode function — we can simply create the truncated exploded dice like this:1

D: {{1..9}:90, {11..20}:9}


Which I'd call close enough:

Anydice will then model up to 8k3, at least fairly close. It will slant slightly lower, and obviously loses out on the extreme highs (which are pretty much 0 anyway).

You can consider the effect for this for pools were anydice is willing to calculate with explodes, say for 5k3:

output [highest 3 of 5dD] named "Truncated 5k3"
output [highest 3 of 5d [explode d10]] named "Exploded 5k3"


Carcer points out you can do the same thing by changeing anydice's explode depth:

set "explode depth" to 1


but I'll stick to the custom die method partly to show it off and because it appears to be slightly faster, but onfortunatly not enough to be give an actual benefit to us here.

• FYI, you can already change the explosion depth by set "explode depth" to 1, and defining a custom die rather than using the built in explosion function doesn't seem to be any more efficient in terms of how large a pool it allows you to consider. – Carcer Mar 20 '20 at 18:06
• Taking this as a base, you get a much smaller margin of error by adding a 2nd. explosion and modeling it like this: anydice.com/program/1a8be - to cope with the larger rolls, you need to remove the last line. – Trish Mar 22 '20 at 15:27
• @Trish: With the default explode depth of 2, [explode d10] is exactly equivalent to {{1..9}:900, {11..19}:90, {21..30}:9}. The only difference to your "truncated" die is that your last range ends at 29 instead of 30. (Effectively, you're conditioning on three consecutive 10s never happening, instead of counting it as 30.) While having 27 instead of 28 possible results for the custom die will surely have some effect on the runtime, It doesn't seem to make any practical difference. (In particular, it seems either way AnyDice can handle XkY for up to X=6 but not for X=7, regardless of Y.) – Ilmari Karonen Mar 23 '20 at 12:59

It looks like AnyDice just can't handle calculating "10k3" even with rerolls limited to just one. So I took my old dice probability calculator written in Python and added a few more features to it.

With the code from this GitHub gist saved as dice_roll.py in the current directory, you can load it into the Python REPL with python -i dice_roll.py and then e.g. calculate and print the distribution of 10k3 (with up to two rerolls by default) in CSV format like this:

exploded_d10 = explode(10, count=2)

for num, prob in sum_roll(exploded_d10, count=10, select=3, ascending=True):
print('%d, %.12g' % (num, 100*prob))


Or just Try It Online!

The results, for various numbers of allowed rerolls per die, look like this:

Looking at the graph, one can see that the first two rerolls make a noticeable difference, but the effect of the later ones is pretty negligible. Which makes sense: for each die, the probability of getting at least n rerolls is 1/10n, so the expected number of dice out of 10 that get 3 or more rerolls is 10/103 = 1/100. And since this expected number of third rerolls is much less than one, it's also approximately equal to the probability of getting even a single third reroll. And the expect number of fourth rerolls, of course, is only 1/1000, and so on.

The Python implementation I wrote handles this problem better than AnyDice for two reasons. The first is simply that it doesn't have the 5 second runtime limit of AnyDice, so (at least when running it locally on your own computer) you can let it run as long as it needs to.

The second reason is that my code is actually a bit smarter than AnyDice, and avoids generating all the possible combinations of the lowest 10 − 3 = 7 dice rolls just to throw them away. This means that, despite Python being a much slower language in general than C++ (which I believe AnyDice is written in), my program actually manages to calculate e.g. 10k3 with max 2 rerolls in only about 0.1 seconds on TIO, significantly faster than AnyDice (which times out).

The code in the gist is actually quite a flexible mini-framework, and can do pretty much anything AnyDice can do (albeit some things more easily than others) and a few things that AnyDice cannot. Some examples:

# basic dice rolls, exploding dice, drop lowest
d6 = make_simple_die(6)                       # d6
sum2d6 = sum_roll(d6, count=2)                # (sum of) 2d6
exp2d6 = explode(sum2d6, count=2)             # [explode 2d6]
output = sum_roll(exp2d6, count=3, select=2)  # [highest 2 of 3d[explode 2d6]]

# statistics (just plain Python, but occasionally useful)
average = sum(n * p for n, p in output)
std_dev = sum((n - average)**2 * p for n, p in output)**0.5

# custom dice are tuples of (value, probability) pairs
dF = tuple((n, 1.0/3) for n in (1, 0, -1))    # fudge die
sum10dF = sum_roll(dF, count=10)              # 10dF

# reverse input die to select lowest instead of highest rolls
rev_d6 = reversed(d6)
lowest = sum_roll(rev_d6, count=4, select=3)  # [lowest 3 of 4d6]

# custom result manipulation example: probability of all dice in 5d6 being equal
yahtzee_prob = 0.0
for roll, prob in dice_roll(d6, count=5):
high = roll[0]  # first element is highest (for normal input dice)
low = roll[-1]  # last element is lowest
if high == low: yahtzee_prob += prob

# dice sides can actually be anything (that can be summed, if using sum_roll)
sqrt_d6 = tuple((n**0.5, p) for n, p in d6)  # sqrt(d6)
sum_sqrt = sum_roll(sqrt_d6, count=3)        # 3d(sqrt(d6))

abcdef = tuple((letter, 1.0/6) for letter in "ABCDEF")
triples = tuple(dice_roll(abcdef, count=3))


The code itself provides basic documentation on how to use the various functions it provides. FWIW, all named arguments in the examples above are optional (with fairly reasonable defaults) and can be either named or given as simple positional arguments, so e.g. sum_roll(d6) and sum_roll(d6, 1, 1) are both equivalent to sum_roll(d6, count=1, select=1).

FWIW, this is getting kind of close to something resembling a reimplementation of AnyDice in Python. I really should consider turning it into a proper Python module with decent documentation.

• Regarding crafting “proper Python module with decent documentation”, if you're open to it, I’d welcome collaboration on my own attempt: posita.github.io/dyce [Deleted and reposted for a link edit. Apologies for the noise.] – posita May 29 at 12:48
• Regarding your implementation of select in dice_roll, I’m having trouble understanding how/why that works or how it’s equivalent to expanding combinations with replacement and selecting from the enumerated results. Do you have (or can you cite) an explanation that might allow me to understand the approach better? – posita May 29 at 15:22
• @posita: Let me try to outline a brief explanation, and you can tell me if it clears up anything. I assume you're familiar with the fact that if the probability of rolling the number $k$ on one die is $p$, then the probability of rolling $k$ exactly $i$ times on $n$ identical dice is $p_i = {n \choose i} p^i (1-p)^{n-i}$? What my algorithm does is take the first side of the die (in the order given in the list, i.e. usually highest or lowest) and let it be $k$ (side in the code) and its probability $p$ (p_side in the code). […] – Ilmari Karonen May 30 at 16:37
• … Then it loops over every possible value of $i$ (keeping in mind that we only care about the $s$ = select highest/lowest rolls, so we only consider the cases where $i<s$ and handle the $i≥s$ case separately later), calculates $p_i$ and then recurses (i.e. calls itself) to generate all possible rolls of the remaining $s-i$ dice (which cannot be $k$, since we already know that $k$ was rolled exactly $i$ times) and their probabilities, which it multiplies by $p_i$ before passing them on to the caller (which might be another instance of itself). – Ilmari Karonen May 30 at 16:45
• BTW, while writing this explanation I realized I could optimize the code a bit: instead of calculating p_select (the probability of rolling $k$ at least $s$ times) as $\sum_{i=s}^n p_i$ at the end of the loop, I can calculate it as $1-\sum_{i=0}^{s-1} p_i$ instead. Which is useful because I already calculate the values of $p_i$ for all $i$ from $0$ to $s-1$ in the previous loop anyway, so I can just sum them up and subtract the result from $1$. – Ilmari Karonen May 30 at 16:57