My friend says you can do this a more theoretical route using an "order statistic" (which I'm not familiar with).

I wrote a program in Scala to calculate this. This should be 100% accurate.

Basic logic:

Iterate every permutation of die values (total of 6^4 permutations)

For each permutation calculate the sum minus lowest. Increment corresponding bucket by one.

Use the permutation counts to calculate expected value.

In the results below, there are 21 permutations for a net value of 18 with a probability of 21/(6^4)=0.016.

Expected value is calculated as normal for discrete variables: (3*1 + 4*4 + 5*10 + 6*21 + ... + 18*21)/(6^4)

My end result is 12.245. This should be 100% accurate.

0: 0: 0.000

1: 0: 0.000

2: 0: 0.000

3: 1: 0.001

4: 4: 0.003

5: 10: 0.008

6: 21: 0.016

7: 38: 0.029

8: 62: 0.048

9: 91: 0.070

10: 122: 0.094

11: 148: 0.114

12: 167: 0.129

13: 172: 0.133

14: 160: 0.123

15: 131: 0.101

16: 94: 0.073

17: 54: 0.042

18: 21: 0.016

expected = 12.245

Code:

object Probability4d6l {
def main(args: Array[String]) {
val permCounts = new Array[Int](19);
val numpermutations = 6 * 6 * 6 * 6;
for (d1 <- 1 until 7) {
for (d2 <- 1 until 7) {
for (d3 <- 1 until 7) {
for (d4 <- 1 until 7) {
val lowest = math.min(d1, math.min(d2, math.min(d3, d4)));
val value = d1 + d2 + d3 + d4 - lowest;
permCounts(value) += 1;
}
}
}
}
var expected: Double = 0.0;
for (i <- 0 until permCounts.length) {
expected += i * permCounts(i).asInstanceOf[Double] / numpermutations;
println("%d: %d: %.3f".format(i, permCounts(i), permCounts(i).asInstanceOf[Double] / numpermutations));
}
println("expected = %.3f".format(expected));
}
}