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 '14 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 '14 at 13:06
• More general version of the question 1. Almost the same question 2. – Affine Aug 5 '14 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 '14 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) – Glen_b Aug 5 '14 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.