probability of subset contained in dice rolls of custom dice with duplicate symbols I play a game where you roll 24 dice, 6 of 4 different colors with each color having a different set of numbers and math operations
these are the dice colors and possible symbols, there are 6 dice of each color, 24 total:
red = 0, 1, 2, 3, +, -
blue = 0, 1, 2, 3, x, /
green = 4, 5, 6, -, x, ^
black = 7, 8, 9, +, /, √
In the game you set a goal of up to 6 dice and the rest of the game involves manipulating the dice to find something equivalent to the goal, but I'm only interested in the probability of setting a goal.
My question:

If I have a subset of symbols that I need say, (x, 1, 1, 4, +, 3), call it $D$. What is the probability that $D$ will exist in each roll of the set of all rolls of 6 red, blue, green, and black dice?
This problem is particularly tricky because some symbols repeat, like x is on blue and green cubes, 1 is on red and blue cubes, + is on red and black cubes.

I've tried to solve this problem myself by brute forcing it in python, but computing 6^24 dice rolls is not feasible for me, and I also tried computing each color as a separate instance, but I couldn't figure out how to find the probability.
I don't know much probability theory, I've taken statistics a while ago, and have looked up numerous questions on this forum; but I either don't know the right question to ask, or I don't find it already asked.
 A: The problem can be solved by using a multivariate generating function, using Mathematica.
As a first step, let's encode the unusual symbols on the dice so they will easier to work with on a computer. 
OP symbol      +    -    x    ^    /    √
Encoded as     10   11   12   13   14   15

With this change, the four colors of dice become
red    0   1   2   3  10  11
blue   0   1   2   3  12  13
green  4   5   6  11  12  14
black  7   8   9  10  13  15

If we use the variable $x_i$ to track the number of occurrences of face $i$, the generating function for a red die is $x_0+x_1+x_2+x_3+x_{10}+x_{11}$.  But notice that we don't really care about faces 0, 2, or 11 since they aren't in $D$, so we might as well lump their associated variables into a single "don't care" variable, $y$.  With that change, the GF for a red die is
$$f_{red} = x_1 + x_3 + x_{10} + 3 y$$
Similarly, the GFs for the other colors of dice are
$$\begin{align}
f_{blue} &= x_1+x_3+x_{12}+3y \\
f_{green} &= x_4+x_{12}+4y \\
f_{black} &= x_{10}+5y
\end{align}$$
The GF for the roll of six dice of each color is
$$g=(f_{red} \cdot f_{blue} \cdot f_{green} \cdot f_{black})^6$$
What this means, for example, is that the number of ways to roll one $1$, two $3$'s, four $12$s and 17 "don't cares" in six rolls of the four dice is the coefficient of $x_1 x_3^2 x_{12}^4 y^{17}$ when $g$ is expanded.  (Note that the exponents must sum to $24$.)
There are $6^{24}$ possible outcomes, all of which we assume are equally likely. Rather than counting the cases with at least two $1$s and at least one each of $3$, $4$, $10$ and $12$, it seems that it might be more efficient to count the cases which do not qualify and subtract from the total.  I.e., we want to sum the coefficients of 
$$x_1^{i_1} x_3^{i_3} x_4^{i_4} x_{10}^{i_{10}} x_{12}^{i_{12}} y^j$$
in $g$ where $i_1 < 2$ or at least one of $i_3, i_4, i_{10}$ or $i_{12}$ are less than $1$, where $j = 24-i_1-i_3-i_4-i_{10}-i_{12}$.  When we do this in Mathematica, we find the sum is 
$$n = 3467290195987632218 \approx 3.46729 \times 10^{18}$$
So the probability asked for in the OP, that $D$ is contained in a set of 6 rolls each of the four colors of dice, is
$$\frac{6^{24}-n}{6^{24}} \approx \boxed{0.268254}$$
A: The probability is $$\frac{635545571166992339}{2369190669160808448} \approx 0.268.$$
You can compute this by dynamic programming.
For a given multiset of symbols $S$ and a set of dice $D$, we will compute the probability $p(S,D)$ that all symbols in $S$ come out in a throw of $D$ as follows:


*

*Base case: $D$ is empty. If $S$ is empty then the probability is 1, otherwise it's 0.

*If $D$ is not empty, choose $d \in D$ arbitrarily. Go over all possible realizations $x$ of $d$. For each of them, if $x$ in $S$, let $q_x = p(S \setminus \{x\}, D \setminus \{d\})$, and otherwise, let $q_x = p(S, D \setminus \{d\})$. The answer is $\sum_x q_x/|d|$ (in your case, $|d| = 6$).


The easiest way to compute this efficiently is using memoization, which is how I computed the number above.
