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
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?