Doug’s MATLAB Video Tutorials

April 10th, 2009

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?
--------------------
References:
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. 

Here is a link to the slides from an earlier talk by Jeff Buzen on the same topic:

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

19:19 UTC | Posted in Topic: Puzzler | Permalink |

You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

12 Responses to “MATLAB Puzzler: Professor and the Rain”

1. Oliver Woodford replied on April 11th, 2009 at 11:02 UTC :
   

   
   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

2. Petter replied on April 11th, 2009 at 11:20 UTC :
   

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

    
   >> puzzler(0.5)
   
   ans =
   
       0.1000
    

3. Oliver Woodford replied on April 11th, 2009 at 11:37 UTC :
   

   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

4. Petter replied on April 11th, 2009 at 11:43 UTC :
   

   There is an analytical solution. It is possible to analytically calculate the left eigenvector of the matrix A in my code.

5. Petter replied on April 11th, 2009 at 11:59 UTC :
   

   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).'
   %Normalize
   Sol = Sol / (sum(Sol))
   
   %Solution
   Sol(1) * p
    

   I get (1-p)/(3-p)*p as the analytical formula.

6. Oliver Woodford replied on April 11th, 2009 at 12:47 UTC :
   

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

7. Tim Davis replied on April 12th, 2009 at 19:43 UTC :
   

   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…

   ;-)

8. dhull replied on April 13th, 2009 at 13:22 UTC :
   

   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! :)

9. João replied on April 13th, 2009 at 22:21 UTC :
   

   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
       end
   
       wetcount(i) = w;
   end
   
   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-.')
   xlabel('p')
    

10. dhull replied on April 14th, 2009 at 13:11 UTC :
    

    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

11. Yi Cao replied on April 17th, 2009 at 13:49 UTC :
    

     
    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;
    end
    wetprob = nW / N;
     

12. Babak Taati replied on May 15th, 2009 at 14:18 UTC :
    

    that’s pretty good Petter.

Leave a Reply

<pre class="code">
a = magic(3);
sum(a)
</pre>

These postings are the author's and don't necessarily represent the opinions of The MathWorks.