I'm wondering if there is a neat way to do or set up a simulation to check my answer. I tried brute force, looking at the number of 666a rolls (21) then the number 665a rolls etc. using multinomial coefficients. I got 258 (out of $6^4$)
Any smart insights would be welcome.
We can indeed use multinomial coefficients to solve this problem. We need a sum of 16, 17 or 18 for the three highest dice. The scenarios we should consider, thus include:
$$6-6-6-(1 \ldots 6)$$ $$6-6-5-(1 \ldots 5)$$ $$6-6-4-(1 \ldots 4)$$ $$6-5-5-(1 \ldots 5)$$
Taking into consideration the number of ways these dice can be rolled, we respectively find:
$$5 {4 \choose 1} + 1 {4 \choose 4} = 21$$ $$4 {4 \choose 2, 1, 1} + 1 {4 \choose 2} = 54$$ $$3 {4 \choose 2, 1, 1} + 1 {4 \choose 2} = 42$$ $$4 {4 \choose 2, 1, 1} + 1 {4 \choose 1} = 52$$
The probability thus equals:
$$\frac{21 + 54 + 42 + 52}{6^4} = \frac{169}{1296} \approx 0.130$$