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. 
 A: 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.
A: Exact answer
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.
