Probability that 2 numbers will occur exactly 4 times. A die is rolled 10 times, what is the probability that 2 numbers will occur exactly 4 times? So I know that the probability of rolling any one number exactly 4 times would be 
$$\frac{10!}{4!6!}\cdot \left(\frac{1}{6}\right)^4 \cdot\left(\frac{5}{6}\right)^6,$$
so would the probability that two numbers occur exactly four times be?
$$\frac{10!}{8!2!}\cdot\left(\frac{1}{6}\right)^4\cdot\left(\frac{1}{6}\right)^4\cdot\left(\frac{4}{6}\right)^2$$
 A: No, you have ignored two effects.  The first is that you seem to be trying to calculate the probability that two specific numbers occur four times each, and the problem asks for any two numbers.  The second is that once one number is given to have occured exactly four times, the odds of the other number also doing so are different that the original odds.
The easy way to do this is to ask how many ways to arrange to get two four-timers are there, and divide by $6^{10}$. So:


*

*There are $\binom{6}{2}$ ways to choose the two numbers.

*There are $\binom{10}{4}$ ways to choose  which of the dice landed on the specific first number.

*There are then $\binom{6}{4}$ ways to choose  which of the dice landed on the specific second number.

*There are then $4^2$ ways to choose the rolls of the remaining two dice.
The net probability is 
$$
\frac{1}{6^{10}} \cdot \frac{6!}{2!4!}\cdot\frac{10!}{6!4!}\cdot \frac{6!}{4!2!}\cdot 4^2
$$
A: After your clarifications, here's an approach to compute favorable ways:
[ Choose faces to come up ] $\times$ [ Permute ]
There are two possibilities:
$4-4-2\;of\;a\;kind: [\binom62\binom41] \times [\frac{10!}{4!4!2!}]$
$4-4-1-1\;of\;a\;kind: [ \binom62\binom42] \times [\frac{10!}{4!4!1!1!}]$
Add up to get favorable ways, and divide by $6^{10}$ to get the final result
A: We can ask about the more  general problem of rolling an $n$-sided die
$m$ times and observing exactly $p$  values that occur $q$ times. This
corresponds to the marked combinatorial species
$$\mathfrak{S}_{=n}(\mathfrak{P}_{=0}(\mathcal{Z})
+ \mathfrak{P}_{=1}(\mathcal{Z}) + \cdots
+ \mathcal{U} \mathfrak{P}_{=q}(\mathcal{Z}) 
+ \mathfrak{P}_{=q+1}(\mathcal{Z}) + \cdots).$$
What we have here is a  sequence of sets containing the positions from
among  the $m$  rolls  where  a value  from  the  $n$ possible  values
occured, where instances of sets of  $q$ positions are marked. We then
obtain the mixed generating function
$$G(z, u) = 
\left(\exp(z) - \frac{z^q}{q!} + u\frac{z^q}{q!}\right)^n.$$
On extracting the coefficient on $[u^p]$ we get
$$m! [z^m] [u^p] G(z, u) =
m! [z^m] {n\choose p} \frac{z^{pq}}{(q!)^p}
\left(\exp(z) - \frac{z^q}{q!}\right)^{n-p}
\\ = m! [z^m] {n\choose p} 
\frac{z^{pq}}{(q!)^p} \frac{z^{(n-p)q}}{(q!)^{n-p}}
\left(\frac{q!\exp(z)}{z^q} -1\right)^{n-p}
\\ = m! [z^m] {n\choose p} 
\frac{z^{nq}}{(q!)^n}
\sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k}
\frac{(q!)^k \exp(kz)}{z^{kq}}
\\ = m! [z^m] {n\choose p} 
\frac{1}{(q!)^n}
\sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k}
(q!)^k \exp(kz) z^{nq-kq}
\\ = m! {n\choose p} 
\frac{1}{(q!)^n}
\sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k}
(q!)^k [z^{m+kq-nq}] \exp(kz) 
\\ = m! {n\choose p} 
\frac{1}{(q!)^n}
\sum_{k=0}^{n-p} {n-p\choose k} (-1)^{n-p-k}
(q!)^k \frac{k^{m+kq-nq}}{(m+kq-nq)!}.$$
Taking into  account that we must  have $m+kq-nq\ge 0$ we  obtain that
$k\ge n -  m/q.$ The problem definition tells us  that $m\ge pq$ (zero
possibilities  otherwise) and  hence  $m/q \ge  p.$  Thus  we  finally
obtain the probability
$$\bbox[5px,border:2px solid #00A000]{
\frac{m!}{n^m} {n\choose p} 
\sum_{k=\lceil n - m/q\rceil}^{n-p} 
{n-p\choose k} \frac{(-1)^{n-p-k}}{(q!)^{n-k}}
\frac{k^{m+kq-nq}}{(m+kq-nq)!}.}$$
Observe that this will yield $${\frac {875}{69984}}$$ for the question
posed by the OP ($n=6, m=10, p=2, q=4$) which incidentally matches the
answer by @MarkFischler.  Here we have used that ${n-p\choose k}$ will
be zero should $\lceil n - m/q\rceil \lt 0$ and $k$ is negative at the
lower end of the range.
There is  some very basic Maple  code to verify the  correctness of
the formula.

PROB :=
proc(n, m, p, q)
    option remember;
    local ind, digits, rolls, res, valset, match;

    if m < p*q then return 0 fi;

    res := 0;

    for ind from n^m to 2*n^m-1 do
        digits := convert(ind, base, n);

        rolls := digits[1..m];
        valset := convert(rolls, `multiset`);

        match := select(ent->ent[2] = q, valset);

        if nops(match) = p then
            res := res + 1;
        fi;
    od;

    res/n^m;
end;

PROB2 :=
proc(n, m, p, q)
    option remember;
    local ind, digits, rolls, res, valset, match;

    if m < p*q then return 0 fi;

    res := 0;

    for ind from n^m to 2*n^m-1 do
        digits := convert(ind, base, n);

        rolls := digits[1..m];
        valset := convert(rolls, `multiset`);

        match :=
        add(`if`(ent[2] = q, 1, 0), ent in valset);

        if match = p then
            res := res + 1;
        fi;
    od;

    res/n^m;
end;


X := (n, m, p, q)->
m!/n^m*binomial(n,p)
*add(binomial(n-p,k)*(-1)^(n-p-k)/(q!)^(n-k)
     *k^(m+k*q-n*q)/(m+k*q-n*q)!, k=ceil(n-m/q)..n-p);

