Wolfram Blog
Jon McLoone

What Shall We Do with the Drunken Sailor? Make Him Walk the Plank!

June 8, 2011 — Jon McLoone, International Business & Strategic Development

Back in 1988 when Mathematica was just a year old and no one in my university had heard of it, I was forced to learn Fortran.

My end-of-term project was this problem: “A drunken sailor returns to his ship via a plank 15 paces long and 7 paces wide. With each step he has an equal chance of stepping forward, left, right, or standing still. What is the probability that he returns safely to his ship?” I wrote a page or so of ugly code, passed the course, and never wrote Fortran again. Today I thought I would revisit the problem.

A drunken sailor returns to his ship

We can code the logic of the sailor’s walk quite easily using separate rules for each case. Firstly, if he is ever on the 16th step or already on the ship, then he is safely on the ship the next time.

step[{x_, 16} | Ship] :=Ship;

If he is ever outside the bounds of the plank, at position 0 or 8, then he falls in the water and stays there:

step[{0, y_} | {8, y_} | Splash] := Splash;

And for all other cases, he steps to one of four new positions.

step[{x_, y_}] := RandomChoice[{{x, y}, {x + 1, y}, {x - 1, y}, {x, y + 1}}];

To simulate an entire journey, we repeatedly step until he reaches one of the named locations.

simulation[start_List] -= Block[{pos = start}, While[Head[pos] === List, post = step[post]]; pos]

We are now able to run a full simulation from a given start point. Let’s be kind to our sailor and have the concierge of the Monte Carlo Casino, where he has been drinking all night, place him in the middle of the end of the plank before leaving him. This is position {4, 1}.

simulation[{4, 1}]

Splash

Using FixedPointList instead of While gives a neat way to visualize his stagger.

Using ListPlot to visualize stagger

Visualization of stagger

Now all we need to do is to run the simulation lots of times and count the successes versus the failures:

Tally[Table[simulation[{4, 1}], {10,000}]]

{{Splash, 8530}, {Ship, 1470}}

This gives us about a 1470(1470+8530) probability of success, or 14.7%, and it took only five lines of code.

Because of the discrete-state nature of the problem, there is another way to look at it that efficiently gives us the probability not only of this journey, but the chances of getting from any one place to any other place in any specific number of steps, and that is to use Markov chains. I didn’t know anything about Markov chains in 1988, so I have no idea how much work it would have been in Fortran, but here is how you might do it in Mathematica.

First we number all of the possible positions. In the 15 × 7 problem there are 105 possible places to stand on the plank. I will designate the water as position 106 and ship as 107 (or more generally water is width × height + 1 and ship as width × height + 2). Now we construct a matrix to represent all the transition probabilities. The element Ma, b contains the probability that we will step from position a to position b in any one step. For example, since the chance of moving from position 1 to position 2 is ¼, M1, 2 = ¼, etc. If we do this right, every row of the transition matrix sums to 1 since there is a 100% probability that we go to some state. Because at any one time we can only step to 4 possible places out of the 107, the resulting matrix is quite sparse, and so using SparseArray will speed up the calculations.

Constructing the matrix is the trickiest part of the problem, but fortunately SparseArray has a rather neat pattern specification that allows us to describe all of the values in a few lines.

Constructing the matrix using SparseArray

Here for example is the transition matrix for a small 3 × 3 plank:

transitionProbabilityMatrix[3, 3] // MatrixForm

Transition matrix for 3 x 3 plank

And here is the transition matrix for our 7 × 15 problem; it has over 10,000 elements, yet all but 422 are zero, making it very efficient as a SparseArray. The nonzero elements follow a similar pattern to the 3 × 3 example.

ArrayPlot[transitionProbabilityMatrix[7, 15]]

Transition matrix for 7 x 15 problem

Here is the first row describing the possible transitions from position 1 (the left-hand edge of the start of the plank). It shows equal chances of going to position 1 (standing still), 2 (stepping right), 8 (stepping forward), and 106 (splash).

transitionProbabilityMatrix[7, 15][[1]] // Normal

Possible transitions from position 1

The important thing about the transition matrix in Markov chains is that Mn represents the transition probability matrix for each start position to end position after exactly n steps. This is my very favorite use of matrices in problem solving. When I learned it, finding powers of matrices was hard work, involving decomposing them into diagonal matrices, but now I can just use MatrixPower. Because MatrixPower supports sparse methods, it is very fast on this problem.

So for example, here are the probabilities for where we will find the sailor who started in the middle of the plank after six steps, showing a 23512 chance (position -2) that he is swimming already.

MatrixPower[transitionProbabilityMatrix[7, 15], 6][[4]] // Normal

Probabilities where we will find the sailor who started in the middle of the plank

This Manipulate lets you visualize his probable position (if still on the plank) along with numeric values of his probability of reaching each final destination.

Visualizing his probable position

So now if we find lim x→∞ Mx then we have the probability that he will make it at some point. I am just going to assume that if he hasn’t made it in 200 steps, then he will probably fall asleep on the plank anyway! So here are the probabilities that he is swimming or home after 200 steps:

MatrixPower[N@transitionProbabilityMatrix[7, 15], 200][[4, -2;;]] // Normal

{0.849989, 0.150011}

In fact, 200 steps is enough to converge to 16 decimal places with this method, making it meaningless to go further in machine arithmetic. Of course in Mathematica, unlike my old Fortran 77 code, we are not limited to machine arithmetic, so here is the exact result for 500 steps:

MatrixPower[transitionProbabilityMatrix[7, 15], 500][4, -2;;]] // Normal

Exact result for 500 steps

This has converged to around 55 decimal places; here are the first 50:

N[%, 50]

The first 50 decimal places

Modeling drunken swimming is a harder problem, which I may revisit when Mathematica can do computational fluid dynamics (CFD).

Leave a Comment

18 Comments


Thales Fernandes

Incredible helpfull! That’s what this blog is about!
Never undestood Markov Chains before this!
Many Thanks.

Posted by Thales Fernandes    June 8, 2011 at 5:58 pm
Nasser M. Abbasi

Thanks for a nice artcle.

But when you say at the end:

“which I may revisit when Mathematica can do computational fluid dynamics (CFD).”

Could you explain what you mean by that?

I mean, why would Mathematica NOT be able to do CFD today? What is missing?

Thanks,
–Nasser

Posted by Nasser M. Abbasi    June 8, 2011 at 10:56 pm
Miroslav Štandera

Very nice post Jon,

especially for students who wonder what the probability is good for :-) I thing it would by nice to add some more visualisation to clarify the result in the first part of your article. I mean to add histograms on the left, right and top of the first diagram to show how many times and where the sailor left the plank. And if the whole diagram is wrapped inside the Manipulate function to show changes along the width and the height it would be an excellent material to study. (Sorry for my bad English.)

Posted by Miroslav Štandera    June 9, 2011 at 2:24 am
Bernard Vuilleumier

Congratulations ! Very interesting simulation. I will use it with my students.

Posted by Bernard Vuilleumier    June 9, 2011 at 4:17 am
Samuel Chen

I wonder if deep down, CFD is nothing but a drunk sailor trying to get home…

Posted by Samuel Chen    June 9, 2011 at 7:50 am
Juan

“Firstly, if he is ever on the 16th step or already on the ship, then he is safely on the ship the next time.” He might be in one of the extremes and step into the water.

Posted by Juan    June 9, 2011 at 11:24 am
William Rummler

@Juan, in this 15×7 example, the 16th step is the ship. I got confused, too, as the phrasing seems to imply “16th step” and “already on the ship” are different states. I think “or already on the ship” should be parenthesized to make this clear.

Posted by William Rummler    June 9, 2011 at 2:02 pm
Hausa

This is fantastic…
You guys at Wolfram are incredible.

Posted by Hausa    June 9, 2011 at 2:25 pm
Mark Lauer

Hi Jon, as always a fascinating post. However an exact answer can be found far more simply and quickly as follows (with no need to worry about convergence):

MatrixPower[
SparseArray[{Band[{1, 1}] -> 3, Band[{2, 1}] -> -1,
Band[{1, 2}] -> -1}, {7, 7}], -15, ConstantArray[1, 7]][[4]]

To briefly explain, since the sailor never moves backward, we can determine all probabilities by working backwards from the ship. The probabilities, x[[1;;7]], of making it to the ship from each position in a given row satisfy the equations 3*x[[i]] == x[[i-1]] + x[[i+1]] + y[[i] where y[[1;;7]] are the corresponding probabilities in the next row. This is because there is a 1/3 chance of heading to each of three adjacent positions on the next non-stationary step. The final row (on board ship) is comprised of probabilities of 1, and each row before that is found by inverting the matrix representing these equations.

Cheers,
Mark

Posted by Mark Lauer    June 10, 2011 at 6:33 am
Jon McLoone

@ Mark
An excellent solution. If I had thought of it, I would have used it. I think there is also a closed form solution to be got from the transition matrix. Something to do with the SVD of the matrix, but I decided the blog was already long enough.

I take comfort from the fact that at least your answer agrees with mine!

Posted by Jon McLoone    June 10, 2011 at 9:31 am
MK

Ah, so this is how to write Mathematica code properly. You can learn a language by reading the documentation, but there’s nothing like seeing an example worked out. Thanks!

Posted by MK    June 10, 2011 at 3:11 pm
Robert

Very nice Jon. Is it possible to get the notebook of this post ? Would be fun to play with it !
Thanks.

Posted by Robert    June 11, 2011 at 1:40 pm
selvarajan

Thank you Now I can flaunt my new found knowledge of Sparse Array among those who still haven’t discovered the power of Mathematica and wonder why would any one buy it.

Posted by selvarajan    June 13, 2011 at 5:44 am
Edson Ferreira

What if we wanted all 10000 paths after 6 steps for example. I mean all sequences like this:

SFRLSF
SFSRRF
RFSRFS
SRFRRR
RFLSRR
RFLSFF

Is it possible a solution for this problem using MatrixPower?

Posted by Edson Ferreira    June 21, 2011 at 2:46 pm
Jon McLoone

@ Edson
The MatrixPower approach only gives you the combined probability of getting from a start to end point over all possible routes. If you want to enumerate the different possible routes, then you might want to do something more like
Permutations[Characters["SSSSSSFFFFFFLLLLLLRRRRRR"], {6}]

Posted by Jon McLoone    June 22, 2011 at 2:56 am
Edson Ferreira

@ John McLoone
What I I wanted is to apply Tally in all 6 steps generated paths in a 10000 paths simulation for instance and see what are the most probable to happen.

Posted by Edson Ferreira    June 22, 2011 at 10:50 am
Anees

Drunken swimming is actually trivial — it’s also called drowning :)
More completely, Lim n->0+ life(n) == NIL.

Posted by Anees    June 27, 2011 at 10:03 pm
Jeff

Very nice integrated view on a simple problem. You should have taught my stochastic processes class!

Posted by Jeff    June 28, 2011 at 8:30 am


Leave a comment

Loading...

Or continue as a guest (your comment will be held for moderation):