Dungeons and Dragons: Ability Score combinations I am trying to figure out how likely it is to roll a certain combination of 6 ability scores in Dungeons and Dragons.
edit:
Dungeons and Dragons is a tabletop role-playing game. You create a character in a fantasy world which has 6 main statistics, Strength, Dexterity, Constitution, Intelligence, Wisdom, and Charisma. Those ability scores determine the powers that a certain character can use and how strong they are.
The ability scores are rolled using 4 6-sided dice (4d6) and then using the 3 highest scores, giving you values $\in [3, 18]$ for the possible values.
My goal is to compute the likelihood of getting a combination like (8, 8, 8, 8, 8, 8) or (15, 14, 13, 12, 10, 8).
/edit
I already worked out the probabilities for rolling a single stat and a question on here confirmed that I did the math right:
3: $\frac{1}{1296}$, 4: $\frac{4}{1296}$, 5: $\frac{10}{1296}$, 6: $\frac{21}{1296}$, 7: $\frac{38}{1296}$, 8: $\frac{62}{1296}$, 9: $\frac{91}{1296}$, 10: $\frac{122}{1296}$, 11: $\frac{148}{1296}$, 12: $\frac{167}{1296}$, 13: $\frac{172}{1296}$, 14: $\frac{160}{1296}$, 15: $\frac{131}{1296}$, 16: $\frac{94}{1296}$, 17: $\frac{54}{1296}$, 18: $\frac{21}{1296}$
Now using the same reasoning I managed to work out the probabilities of combinations based on the sum of ability scores, but only if the outcomes are equally likely. I used Python to compute the probabilities like this:
def probability_of_outcomes(outcomes, repetitions, subsum=slice(None, None, None)):
    probability_dict = {}
    sum_combinations = 0
    for c in combinations_with_replacement(outcomes, repetitions):
        for p in set(permutations(c, repetitions)):
            s = sum(sorted(p)[subsum])
            if s in probability_dict:
                probability_dict[s] += 1
            else:
                probability_dict[s] = 1
            sum_combinations += 1

return probability_dict, sum_combinations


if __name__ == "__main__":
    all_comb, num_combinations_stats = probability_of_outcomes(range(1, 7) 4, subsum=slice(1, None, None))

in the same way I can compute:
abilities_comb, num_combinations_ability = probability_of_outcomes(range(3, 19) 6)

Now with the first computation I can already see that not every number is equally likely, but I am unsure where to add that probability in. Also the approach is not very elegant, as I simply try out every possible combination.
Any ideas?
 A: First order the combination, e.g. $(8, 10, 12, 13, 14, 15)$. Then calculate the exact probability of rolling this combination in this order. This calculation  is simple, just multiply the probabilities of the individual rolls you calculated earlier:
$$p(8) \times p(10) \times p(12) \times p(13) \times p(14) \times p(15)$$
But we're undercounting. After all, it's possible that you rolled $15$ first instead of last. So we multiply by $6!$ to account for all permutations:
$$p(8) \times p(10) \times p(12) \times p(13) \times p(14) \times p(15) \times 6!$$
In this case we're done, because all attributes are unique. When this isn't the case, we're now overcounting. Because in a two-dice example, if you swap the order of rolling two sixes, nothing changes.
So to finish the calculation you need to calculate the number of orderings considering duplicates.
If $a$ is the multiset of attributes, $a_1$ is the first attribute, $a_2$ the second, etc, $|a_1|$ is the number of times the first attribute occurs, and $k$ is the number of distinct attributes, the full formula is:
$$p(a) = p(a_1)^{|a_1|}p(a_2)^{|a_2|}\cdots p(a_k)^{|a_k|} \cdot \frac{k!}{|a_1|!|a_2|!\cdots|a_k|!}$$
Note that in our above example where all attributes were distinct the denominator went away because it was simply $1! \cdot 1! \cdots 1! = 1$.
And in the other extreme, where the multiset only has a single element repeated, like in $(8, 8, 8, 8, 8, 8)$, we get:
$$p(8)^6\cdot \frac{6!}{6!} = p(8)^6$$
As expected.
