by R. Grothmann We play the following game: Throw a dice, and decide, if you want to take the value, or try another throw. If you try again and throw the same or something less, you loose everything. You can repeat that as often as you want. What is the best strategy, and what is the average win? The solution is to work back from the 6 to the 1 and compute the expected win for each throw. The value of a 6 is just 6. For any lower number, we decide to throw again, if the expected value of higher numbers is better than the value we have. We can easily make a function for this.
>function value (n) ...
v=1:n; for i=n-1 to 1 step -1; v[i]=max(v[i],sum(v[i+1:n])/n); end return v endfunction
>value(6)
[ 3.5 3 3 4 5 6 ]
The result is, that we keep 3 to 6, and try again with 1 and 2. The 2 has value 3, since
and the 1 has a better value (!) since
Note that we take the value 3 for the number 2 to compute the value of 1. The expected win is a bit more than 4.
>fracformat; mean(value(6)), longformat;
49/12
What, if we play the same game with a deck of 1000 cards?
>v=value(1000); plot2d(v):
It is clear that v(x)=x for x large enough. The turning point is at 414, where we are better off to try again.
>min(v)
414.595
The expected win is about 627.
>mean(value(1000))
627.093265271
Let us explain this graph. We play the game in the interval [0,1], choosing points there at random. If we get x, and choose to continue, our expected win is the average of v(t) for t>=x.
First of all the turning point s is easy to compute, since v(t)=t for t>=s.
>&integrate(t,t,x,1), &solve(%=x,x), %()
2 1 x - - -- 2 2 [x = - sqrt(2) - 1, x = sqrt(2) - 1] [ -2.41421356237 0.414213562373 ]
The value of s is the second solution. It explains the turning point at 414 we observed in the discrete case for n=1000.
>s=sqrt(2)-1
0.414213562373
Left of s, we can differentiate the equation which defines v and get
with the solutions
We determine lambda, such that the two branches fit.
>lambda=s/exp(-s)
0.626779782174
Now we can make a function of that.
>function v(x) := max(lambda*(exp(-x)),x);
The function looks as expected from the case n=1000.
>plot2d("v",0,1):
The expected win is just lambda as computed above.
>romberg("v",0,s)+romberg("v",s,1)
0.626779782174
We can get the same result in Maxima with symbolic analysis.
>s &= sqrt(2)-1; lambda &= s/exp(-s),
sqrt(2) - 1 (sqrt(2) - 1) E
>&integrate(lambda*exp(-x),x,0,s)+integrate(x,x,s,1) | ratsimp
sqrt(2) - 1 (sqrt(2) - 1) E
It is also possible to compute the discrete results exactly. We seek a discrete function with
It is clear, that v(k)=k for large i.
>&assume(k>0,n>0); sk &= sum(i,i,k+1,n)/n | simpsum | factor
(n - k) (n + k + 1) ------------------- 2 n
We compute the turning value k0. So we search the largest value k, where v(k)>k. In this case, v(k) becomes the sum on the right side of its equation.
>&solve(sk=k,k); k0 &= rhs(%[2])
2 sqrt(8 n + 8 n + 1) - 2 n - 1 ------------------------------ 2
The fraction k0/n has the known limit.
>&limit(k0/n,n,inf)
sqrt(2) - 1
For smaller values than k0, we get
Thus
So it is no surprise, that the exponential function appears. We do in fact get the same result in the limit.
>&(1+1/n)^k0*k0, &limit(%/n,n,inf) | ratsimp, %()
2 sqrt(8 n + 8 n + 1) - 2 n - 1 ------------------------------ 1 2 ((- + 1) n 2 (sqrt(8 n + 8 n + 1) - 2 n - 1))/2 sqrt(2) - 1 (sqrt(2) - 1) E 0.626779782174