combinatorics - Probability that rolling X dice with Y sides and summing the highest Z values is above some value k

Some background: There is an RPG called Legend of the Five Rings, with an interesting dice system. You roll X dice, and keep the highest Z of them. You add those Z dice together. This is phrased as "X keep Z". All of the dice have ten sides, and if you roll a 10, it "explodes" (you roll again). It continues to explode until you roll something other than a ten, adding all the values together, so that one die is worth more than 10.

I'm trying to create my own RPG system, and am looking at different dice options. That said, I'd like to know how to calculate several different probabilities similar to the L5R roll and keep system:

1) What is the probability that X keep Z on Y sided dice is greater than k, with exploding dice?

2) What is the probability that X keep Z on Y sided dice is greater than k, without exploding dice?

3) What is the probability that X keep Z lowest on Y sided dice is greater than k?

Answer

Let $\{S_1, \ldots, S_n\}$ be random vector of i.i.d. outcomes of the die throws. Also let $S_{k \colon n}$ denote the order statistics from that sample. The total score of the $n$-keep-$k$ scheme equals $T = \sum_{i=k+1}^n S_{i \colon n}$.

Because $S_i$ are positive discrete random variables, the technique of finding the probability generating function is the most promising:
$$
\mathcal{P}_T(z) = \mathbb{E}\left(z^T\right)

$$
Once the probability generating function is know, probabilities of possible outcomes can be read off as series coefficients:
$$
\Pr(T=t) = [z^t] \mathcal{P}_T(z)
$$
Moreover, the probabilites $\Pr(T \leqslant t)$ and $\Pr(T > t)$ can be read off as well:
$$
\Pr(T \leqslant t) = \sum_{k=0}^{t} [z^k] \mathcal{P}_T(z) = [z^t] \frac{\mathcal{P}_T(z)}{1-z}
$$
$$

\Pr(T > t) = 1 - \Pr(T \leqslant t) = [z^t] \frac{1 - \mathcal{P}_T(z)}{1-z}
$$

For all of the cases of interest $\mathcal{P}_T(z)$ is a polynomial or rational function in $z$. But finding it requires the knowledge of the joint distribution of the order statistics $\{S_{i\colon:n}\}$.

Using Mathematica, the distribution of the total score can be found as follows:

TotalHighestScoreDistribution[{x_, z_}, dist_] := 
 Block[{v}, 
  TransformedDistribution[Total[Array[v, z]], 

   Distributed[Array[v, z], 
    OrderDistribution[{dist, x}, Range[x - z + 1, x]]]]]

TotalLowestScoreDistribution[{x_, z_}, dist_] := 
 Block[{v}, 
  TransformedDistribution[Total[Array[v, z]], 
   Distributed[Array[v, z], OrderDistribution[{dist, x}, Range[z]]]]]

Here I can provide the answer for explicit choices of the dice systems.

For $6$-keep-$3$ non-exploding 10-sided die:

Similarly, for the keep-lowest scores system:

The exploding die case I did by simulation, mostly because the probability generating function could not be computed in closed form: