Posts Tagged ‘Puzzles’

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?

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 ;-)

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.

Oddly real

March 12, 2011

Consider a complex-valued continuous-time signal, x : \mathbb{R} \to \mathbb{C}, with the property x(t) = - x^*(-t). Show that the even part, x_e, and the odd part, x_o, of signal x are imaginary and real, respectively. Recall that the even and imaginary parts of a signal are defined as follows [1]

x_e (t) := \displaystyle\frac{x(t) + x(-t)}{2}, \qquad{} x_o (t) := \displaystyle\frac{x(t) - x(-t)}{2}.

__________

My solution:

Since x is a complex-valued signal, we can decompose it into its real and imaginary parts

x(t) = x_R (t) + j x_I (t)

where x_R and x_I are real-valued signals. The complex conjugate of x is

x^*(t) = x_R (t) - j x_I (t)

and, from x(t) = - x^*(-t), we know that x(-t) = - x^*(t). The even part is, then

x_e (t) = \displaystyle\frac{x(t) + x(-t)}{2} = \displaystyle\frac{x(t) - x^*(t)}{2} = j x_I (t)

which is purely imaginary, and the odd part is

x_o (t) = \displaystyle\frac{x(t) - x(-t)}{2} = \displaystyle\frac{x(t) + x^*(t)}{2} = x_R (t)

which is real.

__________

References

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

How many votes did Alice get?

January 5, 2011

Last month’s Ponder This challenge was as follows:

In an election, Charles came in last and Bob received 24.8% of the votes. After counting two additional votes, he overtook Bob with 25.1% of the votes. Assuming there were no ties and all the results are rounded to the nearest promille (one-tenth of a percent), how many votes did Alice get?

And here‘s the solution. It was the easiest Ponder This ever.

__________

My solution:

Let a, b, c \in \mathbb{N} be the number of votes gotten by Alice, Bob, and Charles, respectively. If Charles came in last, then a > c and b > c. Since Bob received 24.8% of the votes and he was not the last one, then we deduce that Alice was the first with more than 50% of the votes. Therefore, we have a > b > c.

In the first counting, Bob got 24.8% of the votes. Since the results are rounded off, we have

\displaystyle \frac{b}{a + b + c} = 0.248 + \delta_1

where |\delta_1| < 0.001 is the round-off error. Hence, we obtain

(0.248 + \delta_1) a + (\delta_1 - 0.752) b + (0.248 + \delta_1) c = 0.

Multiplying both sides by 10^3, and letting s_1 := 10^3 \delta_1, we get

(248 + s_1) a + (s_1 - 752) b + (248 + s_1) c = 0.

Note that, since |\delta_1| < 0.001, we have that |s_1| < 1. After counting two additional votes, Charles overtook Bob with 25.1% of the votes. Let c' = c+2 be the new number of votes that Charles got. Hence,

\displaystyle \frac{c+2}{a + b + c + 2} = 0.251 + \delta_2

where |\delta_2| < 0.001 is the round-off error. Thus, we have

(0.251 + \delta_2) a + (0.251 + \delta_2) b + (\delta_2 - 0.749) c = 1.498 - 2 \delta_2

and, multiplying both sides by 10^3, and letting s_2 := 10^3 \delta_2, we get

(251 + s_2) a + (251 + s_2) b + (s_2 - 749) c = 1498 - 2 s_2

where |s_2| < 1. Since there were no ties and these two additional votes were enough for Charles to overtake Bob, then we have c+2 > b > c. As b,c are natural numbers, we conclude that b = c +1.

Hence, we have a linear system of equations in (a,b,c) \in \mathbb{N}^3

\left[\begin{array}{ccc} (248 + s_1) & (s_1 - 752) & (248 + s_1) \\ (251 + s_2) & (251 + s_2) & (s_2 - 749) \\ 0 & 1 & -1 \end{array}\right] \left[\begin{array}{c} a \\ b \\ c \end{array}\right] = \left[\begin{array}{c} 0 \\ 1498 - 2 s_2 \\ 1 \end{array}\right]

where |s_1|, |s_2| < 1. Note that a,b,c must be natural numbers, which means that there should exist (s_1, s_2) \in (-1,1)^2 such that the linear system of equations has a solution in \mathbb{N}^3. Let us relax the system of equations by making s_1, s_2 = 0 and allowing the unknowns to take values in \mathbb{R} rather than \mathbb{N}. Then, we obtain the new linear system of equations

\left[\begin{array}{ccc} 248 & -752 & 248 \\ 251 & 251 & -749 \\ 0 & 1  & -1 \end{array}\right] \left[\begin{array}{c} a \\ b \\ c  \end{array}\right] = \left[\begin{array}{c} 0 \\ 1498 \\ 1  \end{array}\right]

whose solution is (a^*, b^*, c^*) \approx (84.664, 41.168, 40.168). Taking the floor of each component, we obtain a solution in \mathbb{N}^3

(a,b,c) = (84, 41, 40).

Note that 41 / 165 \approx 0.248485, and 42 / 167 \approx 0.251497, meaning that s_1 = 0.485 < 1 and s_2 = 0.497 < 1, which are admissible values. Note also that c+2 > b > c. Hence, all the conditions are satisfied. We thus conclude that Alice got 84 votes.


Follow

Get every new post delivered to your Inbox.

Join 76 other followers