Wolfram Computation Meets Knowledge

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

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

Drunk sailor walk

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

Comments

Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.

18 comments

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

    Reply
  2. 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

    Reply
  3. 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.)

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

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

    Reply
  6. “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.

    Reply
  7. @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.

    Reply
  8. This is fantastic…
    You guys at Wolfram are incredible.

    Reply
  9. 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

    Reply
  10. @ 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!

    Reply
  11. 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!

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

    Reply
  13. 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.

    Reply
  14. 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?

    Reply
  15. @ 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}]

    Reply
  16. @ 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.

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

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

    Reply