Roll 4 Dice, What's the Expected Value of the Sum of the highest 3?

Question

After writing a simulation in python (code at bottom) I realized my calculations are incorrect but can't figure out where I went wrong. Let {$D_1$, $D_2$, $D_3$, $D_4$} be the ordered dice rolls. Let $X = D_2 + D_3 + D_4$.

$$E[X] = E[D_2 + D_3 + D_4] = E[D_2] + E[D_3] + E[D_4]$$

I treated each $D_i$ as an order statistic. The PDF would then be $$f_{D_i}(x)=\frac{4!}{(i-1)!(4-i)!}\frac{1}{6}\left(\frac{x}{6}\right)^{i-1}\left(\frac{6-x}{6}\right)^{n-i}$$ The expectation for each $D_i$ would be, $$E[D_i]=\sum_{x=1}^6x \cdot f_{D_i}(x)$$ With this method I get $D_2 = 2.398$, $D_3 = 3.435$, and $D_4 = 7.022$ for a sum of $12.855$. Approximating it with the code example I get around $\{12.2447, 12.24464, 12.24494\}$ using 100,000 samples which seems like a pretty significant difference.

This question is similar to this one: Finding the expected value of four dice when only the Lowest three count, but I wanted to know why my application of order statistics was wrong.

Code is as follows:

from numpy import array, sort
from numpy.random import uniform

cycles = 100000
total = 0
for i in range(cycles):
    rolls = sort(array(uniform(1,7, size=4), dtype=int))
    total += sum(rolls[1:])

print(total/cycles)

Answer 1

For clarity, allow me to re-introduce the following notation: let $X_i$ be the outcome of the $i$-th roll, $i = 1, 2, 3, 4$. You are interested in determining the pmf (to be more rigorous, it is pmf, NOT pdf) of $X_{(k)}, k = 2, 3, 4$, which should be, for $x = 1, \ldots, 6$:
\begin{align} & P[X_{(k)} = x] = \sum_{j = 0}^{4 - k}\binom{4}{j}\left[\left(1 - \frac{x}{6}\right)^j\left(\frac{x}{6}\right)^{4 - j} - \left(1 - \frac{x}{6} + \frac{1}{6}\right)^j \left(\frac{x}{6} - \frac{1}{6}\right)^{4 - j}\right]. \end{align} This is different from what you presented (which seems to be the pdf of order statistics for continuous r.v.s), which may answer your question "why my application of order statistics was wrong".

---

Alternatively, let me provide a solution which circumvents determining pmf (which is tedious as you may see from above). The goal is to find $E[X_{(2)} + X_{(3)} + X_{(4)}$]. Using $$X_{(2)} + X_{(3)} + X_{(4)} = X_1 + X_2 + X_3 + X_4 - X_{(1)},$$ the calculation could be simplified as you would only need to deal with $X_{(1)}$ (and $E[X_1]$ of course, but that is trivial).

To compute $E[X_{(1)}]$, let's use the formula $E[Y] = \int_0^\infty P[Y > t]dt$ for non-negative r.v. $Y$ (with its discrete version). It follows that \begin{align} & E[X_{(1)}] = \sum_{i = 1}^6P[X_{(1)} \geq i] = \sum_{i = 1}^6P[X_1 \geq i, X_2 \geq i, X_3 \geq i, X_4 \geq i] \\ =& \sum_{i = 1}^6P[X_1 \geq i]^4 = \sum_{i = 1}^6 \left(\frac{6 - i + 1}{6}\right)^4 = 1.755401. \end{align}

On the other hand, \begin{align} E[X_1] = \sum_{i = 1}^6 P[X_1 \geq i] = \sum_{i = 1}^6 \frac{6 - i + 1}{6} = 3.5. \end{align}

Therefore, \begin{align} & E[X_{(2)} + X_{(3)} + X_{(4)}] = E[X_1 + X_2 + X_3 + X_4 - X_{(1)}] \\ =& 4E[X_1] - E[X_{(1)}] = 4 \times 3.5 - 1.755401 = 12.2446, \end{align} agreeing your simulation finding.