I came across a nice answer by Sasha (I believe a Wolfram engineer) on mathematics stack exchange regarding the distribution of the sum of the top $k$ order statistics. (This was an answer to an RPG dice rolling question, where such distributions often arise).
His answer is copied below, edited to just the top $k$ case and for result finalization:
TotalHighestScoreDistribution[{x_, z_}, dist_] :=
Block[{v},
TransformedDistribution[Total[Array[v, z]],
Distributed[Array[v, z],
OrderDistribution[{dist, x}, Range[x - z + 1, x]]]]];
(d = TotalHighestScoreDistribution[{numdice, highest},
DiscreteUniformDistribution[{1, numfaces}]];
pgf = Expectation[z^t, t \[Distributed] d];
cl = CoefficientList[
Series[(1 - pgf)/(1 - z), {z, 0, highest numfaces}], z];
sumdist =
cl // Differences[Reverse@#] & // Reverse //
Append[#, 1 - Tr@#][[highest ;;]] &); // AbsoluteTiming
I found the generating function approach interesting.
The same result can be had more compactly as:
shorter =
PDF[TransformedDistribution[Tr[Array[x, numdice][[-highest ;;]]],
Array[x, numdice][[-highest ;;]] \[Distributed]
OrderDistribution[{DiscreteUniformDistribution[{1, numfaces}],
numdice}, Range[numdice - highest + 1, numdice]]],
Range[highest, highest numfaces]]; // AbsoluteTiming
Depending on the parameters, one or the other is usually markedly faster.
In any case, on larger cases, both flounder. For example, with
{numdice, numfaces, highest} = {12, 6, 6}
They both take on the order of 2 to 4 minutes to complete.
Is there a more efficient means to derive the desired exact distribution using Mathematica?