Probability of ordered dice Suppose you throw five dice, order them from lowest to highest number, then select the third die. What is the probability that the third die was number 2?
Or, in general, when throwing $n$ dice and selecting the $i$th die (sorted from low to high), what is the probability that the number of that die is $k_i$? So the above example was $n=5,i=3,k_i=2$. How do you determine a general formula for finding the probability?
 A: Suppose we have a die with $q$  faces that we roll $n$ times, sort the
outcomes, and ask about the probability  that the value $k$ appears at
position $p.$ We thus have the marked combinatorial class
$$\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}
\textsc{SEQ}_{=k-1}(\textsc{SET}(\mathcal{U}\mathcal{Z})) \times
\textsc{SET}(\mathcal{V}\mathcal{Z}) \times
\textsc{SEQ}_{=q-k}(\textsc{SET}(\mathcal{Z})).$$
This yields the mixed EGF
$$G(z, u, v) = \exp((k-1)uz) \exp(vz) \exp((q-k)z).$$
With $r$ the number of values less than $k$ we obtain from coefficient
extraction
$$\sum_{r=0}^{p-1} [u^r] \exp((k-1)uz) \exp(vz) \exp((q-k)z)
\\ = \sum_{r=0}^{p-1} 
\frac{1}{r!} (k-1)^r z^r \exp(vz) \exp((q-k)z).$$
We must have at least $p-r$ instances of the value $k$ which yields
$$\sum_{r=0}^{p-1} 
\frac{1}{r!} (k-1)^r z^r 
\sum_{j\ge p-r} \frac{z^j}{j!}
\exp((q-k)z).$$
With $n$ rolls of the die we get
$$n! [z^n] \sum_{r=0}^{p-1} 
\frac{1}{r!} (k-1)^r 
\sum_{j\ge p-r} \frac{z^{j+r}}{j!}
\exp((q-k)z)
\\ = n! \sum_{r=0}^{p-1} 
\frac{1}{r!} (k-1)^r 
\sum_{j\ge p-r} \frac{1}{j!}
[z^{n-j-r}] \exp((q-k)z).$$
We must have $n-j-r\ge 0$ or $n-r\ge j$ so we obtain
$$n! \sum_{r=0}^{p-1} 
\frac{1}{r!} (k-1)^r 
\sum_{j=p-r}^{n-r} \frac{1}{j!}
[z^{n-j-r}] \exp((q-k)z)
\\ = n! \sum_{r=0}^{p-1} 
\frac{1}{r!} (k-1)^r 
\sum_{j=p-r}^{n-r} \frac{1}{j!}
\frac{(q-k)^{n-j-r}}{(n-j-r)!}.$$
This yields for the probability
$$\bbox[5px,border:2px solid #00A000]{
\frac{1}{q^n} 
\sum_{r=0}^{p-1} 
{n\choose r} (k-1)^r 
\sum_{j=p-r}^{n-r} {n-r\choose j}
(q-k)^{n-j-r}.}$$
The following  Maple code was used  to verify this formula.  There are
two enumeration routines, one from the problem statement and the other
one incorporating some optimization. The closed formula is implemented
as well and the output matched on all cases that were examined.

with(combinat);

ENUM :=
proc(q, n, p, k)
option remember;
local ind, rolls, res;

    res := 0;

    for ind from q^n to 2*q^n-1 do
        rolls := convert(ind, base, q);
        rolls := sort(rolls[1..n]);

        if rolls[p] = k-1 then
            res := res + 1;
        fi;
    od;

    res/q^n;
end;

ENUM2 :=
proc(q, n, p, k)
option remember;
local res, part, psize, vals, ordpart, rolls;

    res := 0;

    part := firstpart(n);

    while type(part, list) do
        psize := nops(part);

        for vals in choose(q, psize) do
            for ordpart in permute(part) do
                rolls :=
                [seq(vals[blk]$ordpart[blk],
                     blk=1..psize)];

                if rolls[p] = k then
                    res := res +
                    n!/mul(it!, it in ordpart);
                fi;
            od;
        od;

        part := nextpart(part);
    od;

    res/q^n;
end;


X := (q, n, p, k) ->
1/q^n*add(binomial(n, r)*(k-1)^r*
          add(binomial(n-r, j)*(q-k)^(n-j-r),
              j=p-r..n-r), r=0..p-1);

A: We can use the multinomial law to compute this probability. If we split the probability on the number of 2 you draw we have
$$p=\mathbb P(\textrm{exactly one 2 and two 1}) + \mathbb P(\textrm{exactly two 2 and (one or two) 1}) + \mathbb P(\textrm{exactly three 2}) + \mathbb P(\textrm{exactly four 2}) + \mathbb P(\textrm{exactly five 2})$$
$$p=\left(5 \cdot 6 \cdot \left( \dfrac{1}{6} \right)\left( \dfrac{1}{6} \right)^2\left( \dfrac{4}{6} \right)^2 \right) + \left(10 \cdot 3 \cdot \left( \dfrac{1}{6} \right)^2\left(\left( \dfrac{1}{6} \right)^2\left( \dfrac{4}{6} \right)+\left( \dfrac{1}{6} \right)\left( \dfrac{4}{6} \right)^2 \right) \right) + \left(10 \cdot \left( \dfrac{1}{6} \right)^3\left( \dfrac{5}{6} \right)^2 \right) + \left( 5\left( \dfrac{1}{6} \right)^4\left( \dfrac{5}{6} \right)\right) + \left( \left( \dfrac{1}{6} \right)^5\right)$$
$$p=\left( \dfrac{30 \cdot 16}{6^5}\right) + \left( \dfrac{30 \cdot 20}{6^5}\right) + \left( \dfrac{10 \cdot 25}{6^5}\right) + \left( \dfrac{5 \cdot 5}{6^5}\right)+\left( \dfrac{1}{6^5}\right)$$
$$p=\dfrac{1356}{6^5} \approx 0.17438$$
And we get the probability that @Henry found with R. We can also verify it with a simple Python program that gives the same result
from random import randint

ok = 0
nb_draws = 100000
for j in range(nb_draws):
    l = [randint(1,6) for i in range(5)]
    l.sort()
    if l[2] == 2:
        ok += 1

print(ok * 1.0 / nb_draws)

A: A simple general way of solving this is via order statistics. In particular, let  $X \sim \text{DiscreteUniform}(1,6)$ represent our die with pmf $f(x)$:

Given a sample of size $n$ drawn on parent $X$, the pmf of the $i^{\text{th}}$ ordered variable (order statistic) $X_{(i)}$, denoted $g(x)$, is:

where I am using the OrderStat function from mathStatica/Mathematica to automate the nitty-gritties, and where Beta[x,a,b] denotes the incomplete Beta function $\int _0^x t^{a-1} (1-t)^{b-1} dt$. This is the general solution that you seek. For instance, in the specific case provided, with $n = 5$, $i = 3$, and $x = 2$, the exact $P(X_{(3)}=2)$ is:

... which is $\approx$ 0.174383. As disclosure, I should add that I am one of the authors of the software used above.
