Archive for the ‘Puzzles’ Category

Fair play

March 31, 2012

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 p is the probability that the coin shows heads, give the best rational approximation to p with a denominator having fewer than 10 digits.

Here is the solution.

__________

My solution #1:

Let \mathbb{P} (H) = p 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 p < 1/2. The probability that the coin toss will first produce “heads” on the n-th trial (where n \geq 1) is (1-p)^{n-1} p. Let us define \mathcal{I} := \{1, 14, 15, 19, 20, 23\}. The probability that the computer will win is thus

\displaystyle \sum_{n \in \mathcal{I}} (1-p)^{n-1} p.

We now assume that the computer knows the exact value of p 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 1/2. Hence, we obtain the equation

p + (1-p)^{13} p + (1-p)^{14} p + (1-p)^{18} p + (1-p)^{19} p +(1-p)^{22} p = \frac{1}{2}

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

p = 0.499905242928390.

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

p \approx \displaystyle\frac{69693988}{139414397}

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

p =0.49990524292838950452962194522516.

We modify the Haskell script, using the new value of p 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

p \approx \displaystyle\frac{61031369}{122085875}

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 p which I obtained using MATLAB is quite different from the value of p used by the webmaster

p = 0.499905242928506776611191007553.

Which value of p is “more correct”? I do not know. Do you?

Wythoff sequences

August 21, 2011

On November 20, 2007, the Putnam problem of the day was:

The set of all positive integers is the union of two disjoint subsets \{f(1), f(2), f(3), \dots\}, \{g(1), g(2), g(3), \dots\}, where f(1) < f(2) < f(3) < \dots, and g(1) < g(2) < g(3) < \dots, and g(n) = f(f(n))+1 for n=1, 2, 3,\dots. Determine f(240).

Back in November 2007 I could not solve the problem. Now, almost four years later, I will give it another try.

__________

My “solution”:

Let \mathbb{N} := \{1, 2, 3, \dots\} be the set of natural numbers (positive integers). We define F := \{f(n)\}_{n \in \mathbb{N}} and G := \{g(n)\}_{n \in \mathbb{N}}. Note that we must have F \cup G = \mathbb{N} and F \cap G = \emptyset.

Let f(1) = 1. Then, we obtain

g(1) = f(f(1)) + 1 = f(1) + 1 = 2.

Next, let f(2) = 3 and f(3) = 4. Then,

g(2) = f(f(2)) + 1 = f(3) + 1 = 5.

Let f(4) = 6. Then, we obtain

g(3) = f(f(3)) + 1 = f(4) + 1 = 7.

Next, let f(5) = 8 and f(6) = 9. Then,

g(4) = f(f(4)) + 1 = f(6) + 1 = 10.

Next, let f(7) = 11 and f(8) = 12. Then,

g(5) = f(f(5)) + 1 = f(8) + 1 = 13.

Let f(9) = 14. Then, we obtain

g(6) = f(f(6)) + 1 = f(9) + 1 = 15.

Next, let f(10) = 16 and f(11) = 17. Then,

g(7) = f(f(7)) + 1 = f(11) + 1 = 18.

Let f(12) = 19. Then, we obtain

g(8) = f(f(8)) + 1 = f(12) + 1 = 20.

Next, let f(13) = 21 and f(14) = 22. Then,

g(9) = f(f(9)) + 1 = f(14) + 1 = 23.

Let us stop iterating here. So far, we have

\left(f(n)\right)_{n \in \mathbb{N}} = \left(1, 3, 4, 6, 8, 9, 11, 12, 14, 16, 17, 19, 21, 22, \dots\right)

\left(g(n)\right)_{n \in \mathbb{N}} = \left(2, 5, 7, 10, 13, 15, 18, 20, 23, \dots\right).

Consulting the On-Line Encyclopedia of Integer Sequences (OEIS), we find out that \left(f(n)\right)_{n \in \mathbb{N}} and \left(g(n)\right)_{n \in \mathbb{N}} are the lower Wythoff (A000201) and upper Wythoff (A001950) sequences, respectively. These two are Beatty sequences generated by the golden ratio \varphi := \frac{1 + \sqrt{5}}{2}, as follows

f (n) = \lfloor n \varphi \rfloor,     g (n) = \lfloor n \varphi^2 \rfloor

where \lfloor \cdot \rfloor is the floor function. Recall that \varphi^2 = \varphi + 1, which allows us to write g (n) = \lfloor n (\varphi + 1) \rfloor.

The following Python script produces the first twenty terms of the lower and upper Wythoff sequences:

from math import floor
from math import sqrt

# golden ratio
phi = (1 + sqrt(5)) / 2.0

# define function f
def f (n):
    return int(floor(phi * n))

# define function g
def g (n):
    return int(floor((phi+1) * n))

print "Lower and upper Wythoff sequences:\n"
for n in range(1,21):
    print "f(%d) = %d \t g(%d) = %d" % (n, f(n), n, g(n))

This script produces the output:

Lower and upper Wythoff sequences:

f(1) = 1         g(1) = 2
f(2) = 3         g(2) = 5
f(3) = 4         g(3) = 7
f(4) = 6         g(4) = 10
f(5) = 8         g(5) = 13
f(6) = 9         g(6) = 15
f(7) = 11        g(7) = 18
f(8) = 12        g(8) = 20
f(9) = 14        g(9) = 23
f(10) = 16       g(10) = 26
f(11) = 17       g(11) = 28
f(12) = 19       g(12) = 31
f(13) = 21       g(13) = 34
f(14) = 22       g(14) = 36
f(15) = 24       g(15) = 39
f(16) = 25       g(16) = 41
f(17) = 27       g(17) = 44
f(18) = 29       g(18) = 47
f(19) = 30       g(19) = 49
f(20) = 32       g(20) = 52

Finally, we have f (240) = \lfloor 240 \varphi \rfloor = 388.

__________

I am not happy with this “solution”. I obviously cheated when I consulted the OEIS. Supposing one is not acquainted with Beatty sequences, is there any hope of solving the problem? If one happens to be acquainted with Beatty sequences and postulates that

f (n) = \lfloor n r \rfloor,     g (n) = \lfloor n s \rfloor

where \frac{1}{r} + \frac{1}{s} = 1, how does one arrive at the conclusion that r = \varphi and s = \varphi^2? In other words, what would a proper solution look like?

Rainy days

July 7, 2011

A couple of weeks ago, Presh Talwalkar posted this puzzle:

In Mathland, the weather is described either as sunny or rainy. On a sunny day, there is an equal chance it will rain on the following day or be sunny. On a rainy day, however, there is a 70% chance it will rain on the following day versus a 30% chance it will be sunny. How many days is it expected to rain in a 365-day year? Assume the first day is sunny.

__________

My solution:

We will consider a period of N days. Since in our oversimplified world each day is either sunny or rainy, we introduce N random variables W_1, W_2, \dots, W_N taking values in \{S, R\}, where W_k = R means that the k-th day is rainy. We assume that the first day is sunny, i.e., W_1 = S. Hence, \mathbb{P} (W_1 = S) = 1 and \mathbb{P} (W_1 = R) = 0.

Let us now introduce N random variables X_1, X_2, \dots, X_N taking values in \{0, 1\}, where each X_k is defined as follows

X_k =\begin{cases} 1, & W_k = R\\ 0, & W_k = S\end{cases}

Note that the expected value of X_k is \mathbb{E}(X_k) = \mathbb{P} (W_k = R). We introduce a new random variable Y, defined by

Y = X_1 + X_2 + \dots + X_N

where Y is the number of rainy days in a period of N days. We want to find the expected number of rainy days, which is

\mathbb{E} (Y) = \mathbb{E} \left[ \displaystyle\sum_{k=1}^N X_k \right] = \displaystyle \sum_{k=1}^N \mathbb{E} (X_k) = \displaystyle \sum_{k=1}^N \mathbb{P} (W_k = R).

Hence, to compute the desired expected value, \mathbb{E} (Y), we must first compute the probabilities \mathbb{P} (W_k = R) for 1 \leq k \leq N.

Let the probability mass function (p.m.f.) on the k-th day be

\pi_k := \left[\begin{array}{cc} \mathbb{P} (W_k = S) & \mathbb{P} (W_k = R)\end{array}\right]

We know that the p.m.f. on the first day is \pi_1 = \left[\begin{array}{cc} 1 & 0\end{array}\right], as we assume the first day is sunny. Propagating the p.m.f. forwards in time, we obtain \pi_{k+1} = \pi_k P, where

P := \left[\begin{array}{cc} 0.5 & 0.5\\ 0.3 & 0.7\\\end{array}\right]

is the transition matrix. Hence, we have \pi_k = \pi_1 P^{k-1} and, thus,

\mathbb{P} (W_k = R) = \left[\begin{array}{cc} 1 & 0\end{array}\right] \left[\begin{array}{cc} 0.5 & 0.5\\ 0.3 & 0.7\\\end{array}\right]^{k-1} \left[\begin{array}{c} 0 \\ 1\\\end{array}\right]

Matrix P has the eigendecomposition P = Q D Q^{-1} and, hence, P^{k-1} = Q D^{k-1} Q^{-1}. We can finally compute \mathbb{P} (W_k = R) by matrix-vector multiplication, a task so excruciatingly dull that I am most happy to delegate it to my Pythonic minion:

from sympy import *

# create symbolic variable
k = Symbol('k')

# define transition matrix
P = Matrix([[0.5, 0.5], [0.3, 0.7]])

# eigendecomposition of the transition matrix
Q = Matrix([[1/sqrt(2), 5/sqrt(34)], [1/sqrt(2), -3/sqrt(34)]])
D = Matrix([[1, 0], [0, 0.2]])

# check if eigendecomposition is correct
print Q * D * Q.inv() - P

# create vectors
e1 = Matrix(2, 1, [1, 0])
e2 = Matrix(2, 1, [0, 1])

# compute and print probability that the k-th day is rainy
print e1.T * Q * Matrix([[1, 0], [0, 0.2**(k-1)]]) * Q.inv() * e2

This Python / SymPy script gives us

\mathbb{P} (W_k = R) = \frac{5}{8} \left(1 - 0.2^{k-1}\right)

and, finally, we obtain the desired expected value

\mathbb{E} (Y) = \displaystyle \sum_{k=1}^N \mathbb{P} (W_k = R) = \displaystyle \sum_{k=1}^N \frac{5}{8} \left(1 - 0.2^{k-1}\right) = \displaystyle\frac{5}{8}\sum_{k=0}^{N-1} \left(1 - 0.2^k\right)

which is equal to

\mathbb{E} (Y) = \displaystyle\frac{5 N}{8} - \frac{25}{32} (1 - 0.2^N) \approx \displaystyle\frac{5 N}{8} - \frac{25}{32}.

For N = 365, we get \mathbb{E} (Y) \approx 228 rainy days per year. Mathland is certainly not in Southern California ;-)

Dealing cards

April 25, 2011

Imagine that you have a deck of three cards: 2, 3, and 4. You shuffle the deck and deal out the three cards. What is the probability that:

  • the second card is even given that the first card is even?
  • the first two cards are even given that the third card is even?
  • the second card is even given that the first card is odd?
  • the second card is odd given that the first card is odd?

(this is problem 1.5.5 in Yates & Goodman [1])

[ source ]

__________

My solution:

Our experiment consists of shuffling the deck, dealing out the cards, and observing the sequence in which the three cards are dealt. Hence, each outcome or observation is a 3-tuple

(\text{1st card}, \text{2nd card}, \text{3rd card})

and the sample space, the set of all observations, is thus

S = \{(2,3,4), (2,4,3), (3,2,4), (3,4,2), (4,2,3), (4,3,2)\}.

Note that |S| = 3! = 6. We assume that the six observations are equally likely, i.e., \mathbb{P}(s) = 1/6 for all s \in S. Let E_i denote the event that the i-th card dealt is even numbered, and let \mathcal{O}_i denote the event that the i-th card dealt is odd numbered. Recall that an event is a subset of the sample space S and, thus, the event space (the set of all events) is 2^S. Hence, we have

E_1 = \{(2,3,4), (2,4,3), (4,2,3), (4,3,2)\}

E_2 = \{(2,4,3), (3,2,4), (3,4,2), (4,2,3)\}

and E_1 \cap E_2 = \{(2,4,3), (4,2,3)\}. Hence

\mathbb{P}(E_1) = \mathbb{P}\left(\displaystyle\bigcup_{s \in E_1} \{s\}\right) = \displaystyle\sum_{s \in E_1} \mathbb{P}(s) = 4/6

and \mathbb{P}(E_1 \cap E_2) = 2/6. Thus, the probability that the 2nd card is even given that the 1st card is even is the following

\mathbb{P}(E_2 \mid E_1) = \displaystyle\frac{\mathbb{P}(E_1 \cap E_2)}{\mathbb{P}(E_1)} = \displaystyle\frac{2/6}{4/6} = \frac{1}{2}.

Interpretation: if the 1st card is even numbered, after the 1st card is dealt (and observed) there will be only two cards left in the deck, one even numbered and another one odd numbered, and these two are equally likely. Hence, we conclude that the probability that the 2nd card is even given that the 1st card is even is 1/2.

Since there are only two even numbered cards in the deck, we have E_1 \cap E_2 \cap E_3 = \emptyset. Therefore, the probability that the first two cards are even given that the 3rd card is even is

\mathbb{P}(E_1 \cap E_2 \mid E_3) = \displaystyle\frac{\mathbb{P}(E_1 \cap E_2 \cap E_3)}{\mathbb{P}(E_3)} = \displaystyle\frac{\mathbb{P}(\emptyset)}{\mathbb{P}(E_3)} = 0.

Interpretation: if the 3rd card is even numbered, then the two cards previously dealt must be odd and even numbered (not necessarily in that order, of course), as we only have two even numbered cards in the deck. Thus, it’s impossible to get three even cards and \mathbb{P}(E_1 \cap E_2 \mid E_3) = 0.

Note that \mathcal{O}_1 =\{(3,2,4), (3,4,2)\} \subset E_2. Hence, we have that E_2 \cap \mathcal{O}_1 = \mathcal{O}_1. Therefore, the probability that the 2nd card is even given that the 1st card is odd is

\mathbb{P}(E_2 \mid \mathcal{O}_1) = \displaystyle\frac{\mathbb{P}(E_2 \cap \mathcal{O}_1)}{\mathbb{P}(\mathcal{O}_1)} = \displaystyle\frac{\mathbb{P}(\mathcal{O}_1)}{\mathbb{P}(\mathcal{O}_1)} = 1.

Interpretation: if the 1st card is odd numbered, then only even numbered cards remain in the deck. Thus, we know for sure that the 2nd card will be even if the 1st is odd.

Since we only have one odd numbered card in the deck, we can’t have an observation / outcome with two odd numbered cards, i.e., \mathcal{O}_1 \cap \mathcal{O}_2 = \emptyset. Hence, the probability that the 2nd card is odd given that the 1st card is odd is

\mathbb{P}(\mathcal{O}_2 \mid \mathcal{O}_1) = \displaystyle\frac{\mathbb{P}(\mathcal{O}_1 \cap \mathcal{O}_2)}{\mathbb{P}(\mathcal{O}_1)} = \displaystyle\frac{\mathbb{P}(\emptyset)}{\mathbb{P}(\mathcal{O}_1)} = 0.

Interpretation: if the 1st card is odd numbered, then only even numbered cards will remain in the deck. Thus, it’s impossible for the 2nd card to be odd when the 1st one is odd.

__________

References

[1] Roy D. Yates, David J. Goodman, Probability and stochastic processes: a friendly introduction for electrical and computer engineers, 2nd edition, Wiley, 2004.

A convoluted convolution

April 2, 2011

Consider the signal x(t) = 1/t. What is the convolution of x with itself?

__________

My solution:

From the definition of convolution [1], we have

(x * x) (t) = \displaystyle \int_{-\infty}^{+\infty}\frac{d \tau}{\tau (t - \tau)} = \displaystyle \frac{1}{t} \int_{-\infty}^{+\infty} \left(\frac{1}{\tau} + \frac{1}{t - \tau}\right) d \tau.

Unfortunately, the latter integral is plagued by convergence problems. To make matters worse, x(t) is not even defined at t = 0. Fortunately, the former integral is a familiar one: it’s a scaled Hilbert transform of x. It is well-known [1] that the Fourier transform of 1 / \pi t (the impulse response of a Hilbert transformer) is the following

\mathcal{F} \left(\displaystyle\frac{1}{\pi t}\right) = - j \,\text{sgn} (w) =\begin{cases} -j, & \omega > 0\\ +j, & \omega < 0\end{cases}

where \text{sgn}(\cdot) is the signum function. Thus, we have that the Fourier transform of x is

\mathcal{F} (x) =  -j \pi \, \text{sgn} (w) = \begin{cases} -j \pi, & \omega > 0\\ +j \pi, & \omega <  0\end{cases}

Convolution in the time domain corresponds to multiplication in the frequency domain. Hence, the Fourier transform of the convolution is

\mathcal{F} (x * x) = \left(\mathcal{F} (x)\right)^2 = \left(-j \pi \, \text{sgn} (w)\right)^2 = - \pi^2

and, taking the inverse Fourier transform, we finally obtain

(x * x) (t) = \mathcal{F}^{-1} (- \pi^2) = -\pi^2 \delta (t)

where \delta (t) is the Dirac delta. This is a prime example of a problem that Signal Processing professors in engineering departments tend to love, but that drives (some) mathematicians up the wall.

__________

References

[1] Alan V. Oppenheim, Alan S. Willsky, S. Hamid Nawab, Signals & Systems, 2nd edition, Prentice-Hall, 1996.


Follow

Get every new post delivered to your Inbox.

Join 76 other followers