# Probability for completely unfair stats when rolling ability scores using 4d6 drop lowest for D&D

Players of a certain TRPG have characters with 6 ability scores, each ability score ranging from 3-18. One method of generating those is by rolling 4d6 drop lowest. That means four six-faced-dice are rolled, and the three highest results are added.

What's the probability that, given 5 players, one player will have a highest ability score equal or lower than the lowest ability score of another player?

The related question here shows how to get the distribution of 4d3 drop lowest, but how do I get from there to an answer of my question above?

A good answer would explain the result in a way that a statistics novice can follow.

• While this is not a statistical answer, it's worth noting that anydice.com provides a calculator for precisely this purpose.
– Sycorax
Aug 5, 2014 at 13:01
• An answer that contains a script for anydice that shows me the answer to my questions would be accepted (if it was the best or only answer, of course)
– Mala
Aug 5, 2014 at 13:06
• More general version of the question 1. Almost the same question 2. Aug 5, 2014 at 13:10
• @Affine The answers to those are not understandable for a novice. I am asking here because I would like a nice explanation that people who are not otherwise into statistics can follow (and that I can link to). I will update the question to clarify that.
– Mala
Aug 5, 2014 at 13:16
• While these should be possible to do analytically, I'd be inclined to answer via simulation. However, beware of the issue of calculating probabilities for events specified post hoc -- i.e specified after they occur -- as if they'd been specified before the fact. (see the wheelbarrow-full-of-dice example here) Aug 5, 2014 at 13:17

Unlike in answering your other question, I'm just going to do this with Monte Carlo sampling. It wouldn't be super-hard to solve it out explicitly, but it's late. :p If that's what you're looking for let me know and maybe I'll do it out this weekend.

We can get an explicit pmf for $X_{ij}$ by brute force, as before. Now that we have that, let's just roll up 1,000,000 parties real quick, each with 5 players, each with 6 ability scores. (Computers are nice.)

Using scipy and the score_pmf variable from last time:

from scipy import stats
ability = stats.rv_discrete(
name="ability", values=(np.arange(score_pmf.size), score_pmf))
samps = ability.rvs(size=(1000000, 5, 6))
min_max_scores = samps.max(axis=2).min(axis=1)
max_min_scores = samps.min(axis=2).max(axis=1)
print('{:%}'.format(np.mean(max_min_scores > min_max_scores)))


printed out 1.002000%, 0.969900%, 0.974200% for me when I repeated it three times. So "about 1%" seems like a good estimate.

If you want to be a little more rigorous, this corresponds to it happening 10,020 + 9,699 + 9,742 = 29,461 times in 3,000,000 trials, which according to an Agresti-Coull binomial confidence interval means that we can be 99% sure it happens between 0.967% and 0.997% of the time.

This is the easiest way to answer any question about dice rolling like this, though it doesn't provide any statistical insight or generalization.

### Step 1: 4d6 keep highest 3

The problem of keeping the highest can be computed through a variety of methods (example), but brute-force enumerating all $$6^4 = 1296$$ possibilities is trivial for a computer as well. However you prefer to do it, this should be the first step.

import hdroller

one_ability = hdroller.d6.keep_highest(4, 3)


### Step 2: Lowest and highest ability scores for a single player

With the distribution of 4d6kh3 in hand, we then proceed to find the joint distribution of the lowest and highest ability scores for a single player. We can do this inductively: given the distribution after the first $$n$$ scores, we can take the joint distribution with one additional (independent) score and use this to compute the distribution for $$n+1$$ scores. (You could probably still get away with brute-forcing all $$16^6 >$$ 16 million possibilities here, but that's not so interesting.)

def append_score(current, score):
# As we see each score, we track the lowest and highest seen so far.
lo, hi = current
return min(lo, score), max(hi, score)

# Base case: lo = 18 and hi = 3 which will be overwritten in the first iteration.
one_player = hdroller.Die((18, 3))

# Iteratively compute the result.
num_scores = 6
for i in range(num_scores):
# apply() takes care of iterating over joint outcomes and tracking the weights.
one_player = hdroller.apply(append_score, one_player, one_ability)


### Step 3: Highest low score and lowest high score

We can again use an inductive strategy to find the highest of all low scores and the lowest of all high scores. (Brute-forcing looks worse here, with the raw number of possibilities being $$136^5 >$$ 46 billion.)

def append_player(a, b):
# This time, we track the highest of low scores and lowest of high scores.
return max(a[0], b[0]), min(a[1], b[1])

# Base case: hi_of_lo = 3 and lo_of_hi = 18 which will be overwritten in the first iteration.
result = hdroller.Die((3, 18))

# Iteratively compute the result.
num_players = 5
for i in range(num_players):
result = hdroller.apply(append_player, result, one_player)

print(result.sub(lambda hi_of_lo, lo_of_hi: hi_of_lo >= lo_of_hi))

print(result.sub(lambda hi_of_lo, lo_of_hi: hi_of_lo > lo_of_hi))


Chance that the highest low score will be $$\geq$$ the lowest high score:

Outcome Probability
False 96.108983%
True 3.891017%

Chance that the highest low score will be $$>$$ the lowest high score:

Outcome Probability
False 99.026919%
True 0.973081%

which indeed falls under Danica's 99% confidence interval.

### Notes

The provided code uses my hdroller Python package. You can run it in your browser using this JupyterLite notebook.