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.
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
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
function q = puzzler(p) q = p.*(1-p)./(2-p);
% 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
I get (1-p)/(3-p)*p as the analytical formula.
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
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')
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;