Expected outcome for repeated dice rolls with dice fixing Here is another dice roll question.
The rules

*

*You start with $n$ dice, and roll all of them.

*You select one or more dice and fix them, i.e. their value will not change any more.

*You re-roll the other dice.

*Repeat that until all dice are fixed (after at most $n$ rounds).

The final score is the sum of the values of all $n$ dice. Assume you want to maximize the expected score.
The questions

*

*What is the expected score?

*Can the re-rolling strategy be easily phrased?

*Are there situations (possibly for slightly larger $n$) where you would re-roll a 6?

Thoughts
It seems to be counter-intuitive to re-roll a 6, but it would give you one extra roll for all other dice, so maybe it is worth it? Or is there an argument disproving this hypothesis, even without answering the first two questions?
Further reading
I wrote a narrative about this question and the answers on my blog post.
 A: 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“.
A: 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.
A: After coding up some strategies to see if I could get to Joachim Breitner's list of expectations, I've found one that does work for the 9 expectations listed:


*

*Order the $N$ rolled dice in numerical order.
The last will be fixed, of course.

*Starting with a list, $L$, of just the first smallest, and $k=0$, the number of dice to be rerolled:


*

*$m = n(L)$

*If $\sum L < E(k+m)-E(k)$
i.e. whenever the expected gain is greater than the current sum


*

*Set $k=k+m$

*Clear $L$


*Add the next smallest to $L$

*If you've reached the last (max) dice, stop this loop.

*Otherwise, repeat at 2.1.


*The current $L$, including the last, are fixed, and the $k$ rest are rerolled, with expectation $E(k)$.
i.e. $E(N)=\sum L+E(k)$.


I need to update my test code much to try to test for larger $N$...
