I recently asked this question on MSE: The probability of throwing two dice with each result being contained in a set and although I am happy with the answer, I cannot figure out how this generalizes to $n$ or even 3 dice.
To state the case for $n$ dice: If we have the sets $S_1, S_2, \dots, S_n$ of possible outcomes and we throw $n$ indistinguishable dice (with the same number of sides) simultaneously, what is the probability that each die is contained in different sets, that is, what is the probability that the first is contained in set $S_{i_1}$, the second is contained in $S_{i_2}$, the third in $S_{i_3}$, $\dots$ and the $n$-th in $S_{i_n}$ with $i_1 \neq i_2 \neq i_3 \neq \dots \neq i_n$.
An example for $n=2$: When have 2 dice and two sets $A$ and $B$ being $\{ 1, 2, 4\}$ and $\{ 1, 2, 5\}$, respectively. We will end up with the possibilities of $\{ 1,1 \}, \{ 2,2 \}$ which account for a $\frac{1}{36}$ probability each and $\{ 1,2 \}$, $\{ 1,4 \}$, $\{ 1,5 \}$, $\{ 2,4 \}$, $\{ 2,5 \}$, $\{ 4,5 \}$ accounting for a $\frac{2}{36}$ probability each, for a total chance of $\frac{14}{36}=\frac{7}{18}$. As @user2661923 correctly pointed out in the other question, you could also determine this by calculating $\frac{(2\times 3^2) - 2^2}{36} = \frac{14}{36} = \frac{7}{18}$.
Formula for $n=2$: In general the chance for $n=2$ dice with sets $A$ and $B$ is: $$\frac{2! \# A \# B - \#(A\cap B)^2}{6^2}$$.
Formula for $n=3$: This is where I got stuck, I could not obtain the formula for three dice with the sets $A$, $B$ and $C$. However, I came this far: $$\frac{3! \#A \#B \#C - 3\#A\#(B\cap C)^2 - 3\#B\#(A\cap C)^2 - 3\#C\#(A\cap B)^2 \pm \dots}{6^3}$$ I do not know what should be at the dots.
My question: What is the correct formula for $n=3$ (and for larger $n$)?
I also wrote some python code to check some cases and provide support, mostly so that double work is avoided; everthing can be controlled by changing the variable sets (add anoter list to add a die).
from itertools import product
from sympy.utilities.iterables import multiset_permutations
from collections import Counter
# sets = [[1, 2, 3], [3, 4, 5], [4, 5, 6]]
sets = [[1, 2, 4], [1, 2, 5]]
n_dice = len(correct)
counter = Counter()
for roll in product(range(1, max(max(i) for i in sets)+1), repeat=n_dice):
times = 0
for perm in multiset_permutations(roll):
if [perm[i] in sets[i] for i in range(n_dice)] == [True]*n_dice:
times = 1
if times > 0:
counter[tuple(sorted(roll))] += times
count = 0
for k, v in sorted(counter.items()):
count += v
print(f'total: {count}, dice roll: {k}, multiplicity: {v}')
You can optimize the inner loop by framing the problem as finding a perfect matching on a bipartite graph where the left-side vertexes are the dice and the right-side vertexes are the sets, and edge $\left(i, j\right)$ exists iff the number rolled on die $i$ is in set $S_j$. A perfect matching thus corresponds to an assignment of a unique die to each set as desired. This can then be solved using e.g. the Hopcroft-Karp algorithm.
Since you can generate any bipartite graph in this way, the non-triviality of such algorithms suggests that there is no expression for the general case that is both efficient and "closed-form".
You can reduce the outer loop by replacing the full Cartesian product with iterating through the possible multisets, since the dice are indistinguishable. Unfortunately I don't have any further ideas off the top of my head for the outer loop.