Last February’s Ponder This challenge was the following:
To decide who would play first in a human vs. computer game, the human suggested to flip a coin it produced from its pocket, stating “heads you win, tails I win.” The computer replied, “Humans are very devious. I know that your coin is not fair, and that the probability of heads on that coin is slightly less than one half. To make this perfectly fair, flip the coin repeatedly, and I win if the coin first shows heads on toss number 1, 14, 15, 19, 20, or 23. You win if it first shows heads on any other toss.”
If the computer is correct and is the probability that the coin shows heads, give the best rational approximation to
with a denominator having fewer than 10 digits.
Here is the solution.
__________
My solution #1:
Let be the probability that the coin will land with the heads side up. If the coin is biased in favor of the human player, then
. The probability that the coin toss will first produce “heads” on the
-th trial (where
) is
. Let us define
. The probability that the computer will win is thus
.
We now assume that the computer knows the exact value of and, thus, is able to compute its probability of winning. Since the computer is interested in making this game “perfectly fair”, then his probability of winning should equal the human’s probability of winning, i.e., it should equal
. Hence, we obtain the equation
which we can solve using the following Python + SymPy script:
from sympy import *
# create symbolic variable
p = Symbol('p')
# solve polyomial equation
f = p + p*(1-p)**13+p*(1-p)**14+p*(1-p)**18+p*(1-p)**19+p*(1-p)**22
Solutions = solve(f-0.5,p)
# find (purely) real solutions
RealSols = filter(lambda z: im(z)==0, Solutions)
# print solutions
print "\nAll solutions:"
for sol in Solutions:
print sol
# print real solutions
print "\nReal solutions:"
for sol in RealSols:
print sol
This script outputs the following list of solutions:
All solutions: 0.334938061741437 - 1.03291350457235*I 0.334938061741437 + 1.03291350457235*I 1.97217279050061 + 0.143849563303172*I 1.50662666735084 - 0.797180194029411*I 1.92419089868505 + 0.454390584787364*I 0.0780115260218319 + 0.464383156442821*I 0.00197767511752256 - 0.101887310197153*I 1.27721993255011 + 0.951650156864179*I 1.27721993255184 - 0.951650156863701*I 0.247036637335883 - 0.624592516289942*I 0.625456509765409 - 0.845921918004056*I 0.0780115260218319 - 0.464383156442821*I 0.00197767511752256 + 0.101887310197153*I 0.991590970537258 - 0.923668626671437*I 1.79082570892033 - 0.599036471538201*I 1.92419089870934 - 0.454390584808249*I 1.79082570891289 + 0.599036471546724*I 0.499905242928390 0.6254565097654 + 0.845921918004058*I 0.991590970537348 + 0.923668626671824*I 0.247036637335883 + 0.624592516289942*I 1.50662666734992 + 0.797180194030292*I 1.97217279048996 - 0.14384956326771*I Real solutions: 0.499905242928390
Hence, we have that the probability of heads is
.
We can convert this number to a ratio of integers in Haskell using function approxRational in the library Data.Ratio:
import Data.Ratio -- probability of heads p = 0.499905242928390 -- create (infinite) list of epsilons epsilons = map (10**) [-1,-2..] -- create (infinite) list of ratios ratios = map (approxRational p) epsilons -- list of ratios whose denominators have less than 10 digits ratios' = takeWhile (\x -> denominator x < 10^10) ratios
Running the script above in GHCi, we obtain a list of ratios:
*Main Data.Ratio> ratios' [1 % 2,1 % 2,1 % 2,1 % 2,2386 % 4773,2611 % 5223,2636 % 5273, 2638 % 5277,13189 % 26383,44843 % 89703,166183 % 332429, 393036 % 786221,1405961 % 2812455,6243733 % 12489833, 14906352 % 29818355,69693988 % 139414397]
and we finally obtain the following rational approximation
whose denominator has nine digits. According to the webmaster, the rational approximation I presented above is not the best one. Hence, let us now try to obtain a better one.
__________
My solution #2:
Most likely, the reason why I did not obtain the best rational approximation is that the real solution I obtained was not precise enough. Hence, let us use the following MATLAB command to solve the polynomial equation:
solve('p + p*(1-p)^13+p*(1-p)^14+p*(1-p)^18+p*(1-p)^19+p*(1-p)^22 = 0.5')
This command produces the following list of solutions:
ans =
.19776751175225606651411342151624e-2-.10188731019715301414082803578298*i
.19776751175225606651411342151624e-2+.10188731019715301414082803578298*i
.78011526021831885333432975572559e-1-.46438315644282071005250209738747*i
.78011526021831885333432975572559e-1+.46438315644282071005250209738747*i
.24703663733588276238861533742100-.62459251628994218281562354280746*i
.24703663733588276238861533742100+.62459251628994218281562354280746*i
.33493806174143698737516495548579-1.0329135045723528674007841220342*i
.33493806174143698737516495548579+1.0329135045723528674007841220342*i
.49990524292838950452962194522516
.62545650976540488339104798827398-.84592191800405765193777444233924*i
.62545650976540488339104798827398+.84592191800405765193777444233924*i
.99159097053730697830960408610308-.92366862667162608535284662817351*i
.99159097053730697830960408610308+.92366862667162608535284662817351*i
1.2772199325509563050241902871670-.95165015686389448554291677056666*i
1.2772199325509563050241902871670+.95165015686389448554291677056666*i
1.5066266673486985232104204448828-.79718019402868789409470957216463*i
1.5066266673486985232104204448828+.79718019402868789409470957216463*i
1.7908257089204728277307542114304-.59903647154154201043833983995659*i
1.7908257089204728277307542114304+.59903647154154201043833983995659*i
1.9241908987061350990877558129099-.45439058479360516453156947604463*i
1.9241908987061350990877558129099+.45439058479360516453156947604463*i
1.9721727904901564352190617939256-.14384956327517564484908725760533*i
1.9721727904901564352190617939256+.14384956327517564484908725760533*i
We now obtain a much more precise real solution
.
We modify the Haskell script, using the new value of and a much finer list of epsilons (the precisions):
import Data.Ratio -- probability of heads p = 0.49990524292838950452962194522516 -- create (infinite) list of epsilons epsilons = map ((1.005)**) [-1,-2..] -- create (infinite) list of ratios ratios = map (approxRational p) epsilons -- list of ratios whose denominators have less than 10 digits ratios' = takeWhile (\x -> denominator x < 10^10) ratios
Which yields the following rational approximation
whose denominator has nine digits. Although this rational approximation does not match the one in the official solution, I still got my name on the “people who answered correctly” list. Surprisingly, the value of which I obtained using MATLAB is quite different from the value of
used by the webmaster
.
Which value of is “more correct”? I do not know. Do you?
