## Rainy days

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

from sympy import *