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 days. Since in our oversimplified world each day is either sunny or rainy, we introduce
random variables
taking values in
, where
means that the
-th day is rainy. We assume that the first day is sunny, i.e.,
. Hence,
and
.
Let us now introduce random variables
taking values in
, where each
is defined as follows
Note that the expected value of is
. We introduce a new random variable
, defined by
where is the number of rainy days in a period of
days. We want to find the expected number of rainy days, which is
.
Hence, to compute the desired expected value, , we must first compute the probabilities
for
.
Let the probability mass function (p.m.f.) on the -th day be
We know that the p.m.f. on the first day is , as we assume the first day is sunny. Propagating the p.m.f. forwards in time, we obtain
, where
is the transition matrix. Hence, we have and, thus,
Matrix has the eigendecomposition
and, hence,
. We can finally compute
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
and, finally, we obtain the desired expected value
which is equal to
.
For , we get
rainy days per year. Mathland is certainly not in Southern California ;-)
Tags: Markov Chains, Probability Theory, Puzzles, SymPy
July 8, 2011 at 00:08 |
This very short script computes the desired expected value: