Probability of ordered dice

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);

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)