Stuart’s MATLAB Videos

Watch and Learn

MATLAB Puzzler: Professor and the Rain 12

Posted by Doug Hull,

Being a company full of geeks, a fair amount of brainteasers get passed around. One that recently came around got my attention. I knew there was probably math that governed it elegantly, but I did not know how to do that. I do know how to brute force something in MATLAB, so I wrote a monte carlo simulation.

I am against premature optimization. I think you should write clear, easy to maintain code and then and only then optimize if it really needs it. In this case, I needed to optimize to get the code to run in a reasonable time. Once you all have had a chance to solve this puzzler in the comments, I will make a video showing the progression from code that was written quickly to code that was optimized for speed. Until then, try this puzzler out:

A professor has an umbrella in the office and another one at home.

If it is raining when the professor goes out, he takes an umbrella with him, 
if one is available.  If it is not raining, he doesn't carry any umbrellas 
with him.  Whenever the professor walks with an umbrella, he stays dry. 
Whenever he walks without an umbrella in the rain, he gets wet.

Given p = probability of rain (same on every trip), what is the (steady 
state) probability of getting wet?

I.e., what is the probability of rain & no available umbrella?
I heard about this problem during Jeff Buzen's talk at the joint IEEE/ACM 
meeting, The Improbable Success of Probabilistic Models.
The problem is taken from the Introduction to Probability textbook by MIT 
professors D.P. Bertsekas and J. N. Tsitsklis. 

As always with these puzzlers, be sure to use the proper tags in the comments.

<pre> <code>
[deltaX, deltaY] = puzzler(inMatrix)
%% All the code so someone can just copy and paste it from the comments.
</code> </pre>

12 CommentsOldest to Newest

Oliver Woodford replied on : 1 of 12

function q = puzzler(p)
q = p.*(1-p)./(2-p);
I worked that formula out by hand. I'd be interested to know if the brute force methods agree with it, or whether anybody comes up with a different analytic solution. Oliver
Petter replied on : 2 of 12
Markov-chain solution:
% States
% 1 -- Professor with no umbrellas
% 2 -- Professor with 1 umbrella
% 3 -- Professor with 2 umbrellas
% p -- probability of rain
function prob = puzzler( p )
    A = [0 0 1;  %If the professor has no umbrella, he will for
                 %sure end up with 2 at his destination
         0 1-p p; %If has has 1 umbrella, he will bring it 
                  %with probability p
         1-p p 0]; %If he has 2 umbrellas, he will bring one
                   %of them with prob. p
    dist = [0 1 0]; %prior
    dist = dist * A^1000; %steady-state
    %The probability of getting wet is the probability
    %of rain times the probability of having no umbrellas
    prob = dist(1) * p;
>> puzzler(0.5)

ans =

Oliver Woodford replied on : 3 of 12
Blast - I made a mistake in my analytical solution. I still think it is possible to evaluate, though. I may come back with a correction. Nice solution, Petter. Oliver
Petter replied on : 4 of 12
There is an analytical solution. It is possible to analytically calculate the left eigenvector of the matrix A in my code.
Petter replied on : 5 of 12
Maybe something like this, using the matrix from before:
syms p
A = [0 0 1;  %If the professor has no umbrella, he will for
             %sure end up with 2 at his destination
     0 1-p p; %If has has 1 umbrella, he will bring it 
              %with probability p
    1-p p 0]; %If he has 2 umbrellas, he will bring one
              %of them with prob. p
%Calculate eigenvectors of the transpose
[V E] = eig(A.')

%Find eigenvector with eigenval. 1
ind = find( diag(E) == 1 )
Sol =  V(:,ind).'
Sol = Sol / (sum(Sol))

Sol(1) * p
I get (1-p)/(3-p)*p as the analytical formula.
Oliver Woodford replied on : 6 of 12
OK, I got there in the end and also get p(1-p)/(3-p). The way I did it was write an expression for q(0) as a function of p and q(2), and an expression for q(2) as a function of p and q(0), where q(x) is the probability of having x umbrellas (these expressions are both geometric series). I then substituted out q(2), giving an expression for q(0), which I then just multiply by p to give the answer. I wouldn't call it elegant :)
Tim Davis replied on : 7 of 12
What's the rest of the story? Don't leave us hanging here, with just the bare facts. Is the professor upset at getting wet, and if so does he walk faster? Does he wear jeans or polyester; the latter taking longer to dry and feels worse when he gets wet? Is he bald or not? (Since it's well know that a full head of hair increases the chance of rain, due to an increase static electricity discharged into the atmoshere). Does the professor meet her true love at a rainy bus stop when a passerby offers her his own umbrella? If so, what is the probability that they will have four children, two boys and two girls, as compared with a different gender mix of four children? (Hint: assuming 4 children, the answer is not 50%). Have to inked a contract for the book rights, and when will the DVD come out? Inquiring MATLAB'ers Want to Know... ;-)
dhull replied on : 8 of 12
Everyone, These are much better answers than I came up with. I knew there was cleaner math than what I used. Twenty minutes programming will save you two minutes in a library... :) Later this week, I will post the Monte Carlo code optimization and a graph of the solution. -Thanks, Doug PS, Tim- Pretty sure the professor wore polyester, most of mine did at least! :)
João replied on : 9 of 12
Here is my Monte Carlo like simulation.
function wetprob = puzzler

p = linspace(0.1,0.9,20); % prob. of rain
N = 1e6; % nr. of trips

wetcount = zeros(size(p)); % nr. of times the professor gets wet

for i = 1:numel(p)

    % initial state
    s = 1; % one umbrella available (can be 0, 1 or 2)
    w = 0; % nr. of times the professor gets wet

    q = p(i);
    for n = 1:N
        r = rand < q; % rain?
        w = w + r*(0.5*(s-3)*s + 1);
        s = (1-r)*abs(s-2) + r*(0.5*(1-s)*s + 2);
        %%% equivalent switch ... case structure
        % switch r
        %     case 0
        %         switch s
        %             case 0, s = 2;
        %             %case 1, s = 1;
        %             case 2, s = 0;
        %         end
        %     case 1
        %         switch s
        %             case 0, s = 2; w = w + 1;
        %             case 1, s = 2;
        %             case 2, s = 1;
        %         end
        % end
    wetcount(i) = w;

wetprob = wetcount/N; % simulated probability of getting wet
teoprob = (1-p).*p./(3-p); % theoretical probability of getting wet

plot(p, wetprob, 'b--o', p, teoprob, 'r-.')
dhull replied on : 10 of 12
João, Wow, that is lightning quick. Much faster than what I came up with (video to be posted this week). I think this is the perfect example of optimizing for speed. You left your original code as reference and commentary and managed to come up with an obfuscated but equivalent and faster set of code. In this case, the speed was worth the extra time and effort it takes to maintain the code. Thanks for a great solution! Doug
Yi Cao replied on : 11 of 12
function wetprob=puzzler(P,N)

% P - probability of rain
% N - number of simulated steps
% wetprob - Probability to get wet

s = 1;  % state represents number of umbrellas (0~2). Initially, one in office and one at home
R = rand(N,1)<P;
nW = 0; % number of wet times
for k=1:N
    r = R(k);
    w = r & ~s;         % if rain and no umbrella then wet
    s = 2 - s + r - w;  % next step, 2-s initial, r-w to carry 
    nW = nW + w;
wetprob = nW / N;