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:
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
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.
% 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
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.
There is an analytical solution. It is possible to analytically calculate the left eigenvector of the matrix A in my code.
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.
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 :)
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…
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.
PS, Tim- Pretty sure the professor wore polyester, most of mine did at least! :)
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')
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!
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;
that’s pretty good Petter.