my question is regarding exploding dice. Exploding dice are dice that can be rolled again if the maximum number has been rolled. At the end, the rolls are added up.
For example, for a 6-sided die like this: Roll 1: 6 => roll again. Roll 2: 6 => roll again. Roll 3: 3. Final result: 6 + 6 + 3 = 15
What I want is a function P(X, Y, Z), that returns the probability, that the sum of two dice exceeds a certain threshold, where X = sides of die 1, Y = sides of die 2 and Z = threshold that needs to be exceeded.
For example: The probability, that the sum of the rolls of two exploding 6-sided dice is at least 8 is ~38%. So P(6, 6, 8) = 38% (according to my simulation).
Thank you for helping me! :)
So, let's look at a single die. Without rerolls, this would be a simple uniform distribution, in which every number has the same probability. let us define $d$ as the number of sides of our die, and $t$ as our target number. If we didn't have rerolls as part of our picture, the Probability Mass Function for our die roll would just be:
$$P(X=t)=\frac{1}{d}$$
The fact that we have rerolls means that, if we want to reach numbers higher than $d$, then we need to roll again, and this will multiply the results for the results from $d$ to $2d$ by $\frac{1}{d}$. So if we take rerolls into account the denominator will be $d^{\left \lceil \frac{t}{d} \right \rceil}$ (ie,$d$ to the power of $\frac{t}{d}$ rounded up to the next integer). Thus the PMF for a single exploding dice is:
$$P(X=t)=\frac{1}{d^{\left \lceil \frac{t}{d} \right \rceil}}$$
And this works, because your scenario notes that you can choose to reroll - thus multiples of 6 are still possible. We don't need to worry about summing up a list of probabilities here, because with 1 rerolling die, there's only ever one combination of rerolls that will land you at any given number.
Two dice is mostly straightforward. The probability for any given roll is going to be the probability of the first die's result multiplied by the probability of the second die, which looks like:
$$P(X=t_1, Y=t_2)=\frac{1}{d_1^{\left \lceil \frac{t_1}{d_1} \right \rceil}d_2^{\left \lceil \frac{t_2}{d_2} \right \rceil}}$$
Of course, there are multiple ways two dice can add up to $t$, so we'll need to add the probabilities of all the matching rolls together. If we set $t_1=k$, then $t=t_2+k$, which we can rearrange to $t_2=t-k$. We now have our two die's results in terms of $t$ and $k$, so we can set up a summation to go through all possible combinations by setting our iterator $k=1$, and having it stop when $k=t-1$:
$$P(X+Y=t)=\sum_{k=1}^{t-1}\frac{1}{d_1^{\left \lceil \frac{k}{d_1} \right \rceil}d_2^{\left \lceil \frac{t-k}{d_2} \right \rceil}}$$
And thus, we have our Probability Mass function for an arbitrary $d_1$,$d_2$ and $t$. To have the function calculate the probability of equalling or exceeding $t$, because $t$ has no upper bound it's easier to calculate the probability of not reaching $t$ and then subtracting 1 from that figure. All we have to do is add a second iterator to our function to count up from 2 (our minimum possible $t$) to $t-1$:
$$P(X+Y\geq t)=1-\sum_{i=2}^{t-1}\sum_{k=1}^{i-1}\frac{1}{d_1^{\left \lceil \frac{k}{d_1} \right \rceil}d_2^{\left \lceil \frac{t-k}{d_2} \right \rceil}}$$
Let's test our function on your test case:
\begin{align*} P(X+Y\geq 8)&=1-\sum_{i=2}^{8-1}\sum_{k=1}^{i-1}\frac{1}{6^{\left \lceil \frac{k}{6} \right \rceil}6^{\left \lceil \frac{8-k}{6} \right \rceil}}\\ &=1-\left ( \frac{1}{36}+ \frac{1}{18}+ \frac{1}{12}+ \frac{1}{9}+ \frac{5}{36}+ \frac{1}{6} \right )\\ &=1-\frac{7}{12}\\ &=\frac{5}{12}=0.41\overline{6} \end{align*}
So, around 41.6% - about 3 percentage points off your simulation.