Expected outcome for repeated dice rolls with dice fixing

For standard six-sided dice, $v_{199}-v_{198}\approx6+9\cdot10^{-21}$. Here's the code I used to compute this value. The numbers are represented exactly as rationals, so rounding errors are excluded.

The basic idea is to keep track of two values for each $n$: the value $v_n$ of the game with $n$ dice and the value $w_n$ of the game when preceded by a free roll without $6$s and without the requirement to fix at least one die. For $n\ge6$, we have $v_n\gt5n$, so the free roll becomes worthless and $w_n=v_n$. Thus we only have to calculate $w_n$ and $v_n$ separately up to $n=5$, and as long as the optimal strategy fixes all $6$s, for $n\gt6$ we can calculate $w_n=v_n$ as

\begin{align} w_n &=6^{-n}\left(5^n(w_{n-1}+5)-\sum_{d=1}^4d^n+\sum_{k=0}^{n-1}\binom nk5^k\left(w_k+6(n-k)\right)\right)\\ &=C_n+6^{-n}\left(\sum_{k=0}^n\binom nk5^kw_k-5^n(w_n-w_{n-1})\right)\;, \end{align}

where in the first equation the first term corresponds to the case without $6$s, the second term corrects for the fact that the highest die in that case need not be a $5$, and the sum in the third term is over the number $k$ of non-$6$s; and where

$$C_n=n+6^{-n}\left(5\cdot5^n-\sum_{d=1}^4d^n\right)\;,$$

as defined in Milo Brandt's answer, is the expected value of the dice to be fixed either as $6$s or as the highest die. (This is essentially the modification of Milo's approach to account for the exceptions at small $n$, as suggested in his answer.)

Since the calculation yields $v_{199}-v_{198}\gt6$, the optimal strategy for $200$ dice does not fix the second of exactly two $6$s. Beyond that, the calculation is no longer valid, since it assumes that $6$s are always fixed.

I added this question here as an answer to Examples of apparent patterns that eventually fail.

I'm leaving the original answer up in case the process that led to the solution is of any interest.

Here are some results vor $v_n-v_{n-1}$ for other dice. They don't settle the concrete question for standard six-sided dice, but they indicate why it's so difficult to settle, and they disprove the intuition discussed in the question that the highest number should never be re-rolled.

For three-sided dice with numbers $0,1,1$:

1 : 0.666666666666666
2 : 1.037037037037037
3 : 1.050754458161865
4 : 1.036156412470998
5 : 1.023458163050328
6 : 1.01495109027056
7 : 1.00954958850473
8 : 1.0061461980586
9 : 1.003988393919573

and so on, approaching $1$ from above until at least $n=50$. So in this case $1$s should always be rerolled, except if all dice are $1$s. But yet more interesting, for three-sided dice with numbers $1,2,3$:

 1 : 2.0
 2 : 2.555555555555555
 3 : 2.818930041152263
 4 : 2.925621094345374
 5 : 2.966531248686746
 6 : 2.983754306121206
 7 : 2.991943835462998
 8 : 2.996014709747178
 9 : 2.998033031889117
10 : 2.999026265538337
11 : 2.999514808948095
12 : 2.999756553231309
13 : 2.999877305913408
14 : 2.999938178399316
15 : 2.999969029284208
16 : 2.99998464607431
17 : 2.999992480293932
18 : 2.999996346687621
19 : 2.999998215606247
20 : 2.999999102436438
21 : 2.999999522113371
22 : 2.999999727722627
23 : 2.999999837836287
24 : 2.999999904492389
25 : 2.999999948801448
26 : 2.999999978835654
27 : 2.999999997994604
28 : 3.000000008463403
29 : 3.000000012352717
30 : 3.000000011876478
31 : 3.000000009460664
32 : 3.000000006368981
33 : 3.000000003409548
34 : 3.000000000987628
35 : 2.999999999170625
36 : 2.999999998255086
37 : 2.999999998245222
38 : 2.999999998631135
39 : 2.999999999186528
40 : 2.999999999746
41 : 3.000000000210072
42 : 3.000000000535028
43 : 3.000000000716657
44 : 3.000000000774869
45 : 3.000000000741289
46 : 3.000000000657582
47 : 3.000000000549125
48 : 3.000000000434833
49 : 3.00000000036221
50 : 3.000000000285031

(Scroll down to see the results up to $n=50$.) The values are below $3$ for $1\le n\le27$, then above $3$ for $28\le n\le34$, then below $3$ for $35\le n\le40$ and then above $3$ until at least $n=72$, approaching it from above after a local maximum at $n=44$.

These results were obtained with the same code that reproduces Joachim's and Byron's results for standard dice up to $n=7$, so a bug is unlikely.

The results show that general strategy stealing arguments (as suggested by mercio in a comment; I'd also tried to find one) will not work and any proof will have to be based on the specific numbers on the six-sided dice. (A six-sided die with one low and five high numbers behaves much like the three-sided die with one low and two high numbers.)

P.S.: Here are the results for two-sided dice (a.k.a. coins) with numbers $0,1$:

1 : 0.5
2 : 0.875
3 : 1.0
4 : 1.0078125
5 : 1.00634765625
6 : 1.0043106079101562
7 : 1.0027518272399902
8 : 1.0016952995210886
9 : 1.0010172286929446

and so on, approaching $1$ from above roughly exponentially until at least $n=60$.

I'm now running a longer calculation with standard four-sided dice to see whether they eventually follow the pattern of the three-sided dice. So far, up to $n=50$, the differences seem to be approaching $4$ from below roughly exponentially. But note that for the three-sided dice with numbers $1,2,3$, the differences appear to approach $3$ from below roughly exponentially all the way right up to $n=27$, where the difference from $3$ is only $2\cdot10^{-9}$, and then they decide to exceed $3$ after all; so all bets are off.

$Update$: For standard four-sided dice, $v_{70}-v_{69}\approx4-1.5\cdot10^{-13}$ and $v_{71}-v_{70}\approx4+1.7\cdot10^{-13}$. This suggests that we may eventually have $v_n-v_{n-1}\gt6$ for standard six-sided dice after all; and since the switch occurs at $n=28$ for three-sided dice and at $n=70$ for four-sided dice, it may occur at quite high $n$ well over $100$ for six-sided dice. Perhaps the effort we've been putting into proving $v_n-v_{n-1}\lt6$ should better be directed at making larger values of $n$ numerically accessible.

Not a full answer, but at least some evidence. I calculated the expected value for the few cases:

1:  3.50000 (+3.50000)
2:  8.23611 (+4.73611)
3: 13.42490 (+5.18879)
4: 18.84364 (+5.41874)
5: 24.43605 (+5.59241)
6: 30.15198 (+5.71592)
7: 35.95216 (+5.80018)
8: 41.80969 (+5.85753)
9: 47.70676 (+5.89707)

So for example when you have rolled a 6-5-1-1, it is better to re-roll three dice rather than keep the 5, as the expected value for three is more than 5 larger than for two dice.

The code is this Haskell code. It uses dynamic programming, but for each number of dice goes through all possibilities, hence I stopped at 9 dice:

import Numeric.Probability.Example.Dice
import Numeric.Probability.Distribution (expected)
import Control.Monad
import Data.List
import Text.Printf

probs = map prob [0..]

prob 0 = 0
prob n = expected $ do
    dice <- dice n
    let sorted = reverse $ sort dice
    return $ maximum 
        [ fromIntegral (sum (take m sorted)) +  (probs !! (n - m)) | m <- [1..n] ]

main :: IO ()
main = forM_ (zip3 [1..9] (tail probs) probs) $ \(n, e, p) ->
    printf "%d: %8.5f (+%7.5f)\n" (n::Int) (realToFrac e::Double) (realToFrac (e - p)::Double)

It seems that the differenced are approaching 6 from below. If that is the case, and they never surpass 6, then the answer to the third question is „no“.

This is not an answer, but does something to the problem. I'm not sure it's useful, but maybe someone can find a way to solve/use the identity I end with.

One thought to attack this is to assume that the best strategy is to fix all $6$'s that come up and, if none come up, to fix the highest number that does come up. This is obviously not the best strategy for small $n$ (or when all but a few dice com up as $6$), but it would not be too hard to modify the below machinery to accommodate this change. Basically, then the idea would be to try to solve the modified system and then check if the difference of the expected values ever exceeds $6$ - if they do, this strategy wasn't optimal after all, if not, then it was optimal.

What strikes me as a natural way to attack this is to try to use generating functions. In particular, suppose that $a_n$ is the sequence of expected values of the strategy. We find the recurrence relation: $$a_n = \left[\sum_{i=1}^{n}{n\choose i}(5/6)^{n-i}(1/6)^ia_{n-i}\right]+(5/6)^na_{n-1}+C_n$$ where $C_n$ is the expected value of the dice fixed in that roll. To find $C_n$, note that the number of $6$'s that we expect to roll $n/6$. The probability that the highest number is a $5$ is $(5/6)^n(1-(4/5)^n)$. The probability that it is a $4$ is $(4/6)^n(1-(3/4)^n)$ and so on. If we sum everything, we find that we expect to fix a value of $$C_n = n + 5(5/6)^n - (4/6)^n - (3/6)^n - (2/6)^n - (1/6)^n$$ on a roll with $n$ dice.

Now, we just encode everything in an exponential generating function. That is, we try to reason about the function $$A(x)=\sum_{n=0}^{\infty}\frac{a_nx^n}{n!}.$$ We use exponential generating functions since summing over a binomial coefficient is a natural thing to do in that context, corresponding to multiplying generating functions together. Note that this basically tells us that the $n^{th}$ derivative of our generating function should be $a_n$ at $0$.

For convenience, we need to encode an exponential generating function for $C$ as well, which happens to be $$C(x)=xe^x + 5e^{5x/6} - e^{4x/6} - e^{3x/6} - e^{2x/6} - e^{x/6}.$$

Then, encoding the recurrence for $a_n$ as a generating function gives $$A(x)=(e^{x/6}-1)A(5x/6)+\int_{0}^{5x/6}A(t)\,dt + C(x).$$ I don't see much hope in actually solving this - it's just too ugly. Moreover, to correct the expected values for the first $k$ terms, which correspond to a different strategy one would have to add a polynomial of degree $k$ to the right hand side. One would have to replace the instance of $A(5x/6)$ with that plus some other polynomial to deal with the fact that the strategy is different when all but a few dice are $6$'s. Neither of these really makes things worse - they just amount to adding some particular function to the right hand side - but they don't make it easier either.

You can check that this expresses what we desire by differentiating it at $x=0$ however many times you want. You'll find that each derivative of $A$ at $0$ will be written in terms of derivatives of lower orders, so at least it does what it's supposed to.

The bright side is that we can work with actual functions, since, due to the bound $a_n\leq 6n$, we have that the sum for $A$ converges to an entire function $\mathbb C\rightarrow\mathbb C$. So, it's possible that some clever application of complex analysis using the above identity does something useful, though I can't think of any criterion that would tell us whether $6$ actually is a bound for the difference of two derivatives or not.