UC SCI 111, v2.0
Notebook 7 - Markov Chains

Initializations

Needs["Graphics`Graphics`"]

Off[General::spell1]

histogram[dat_, opts___Rule] := 
      Histogram[dat, opts, 
        HistogramCategories -> Table[i -  .5, {i, Min[dat], Max[dat] + 1}], 
        Ticks -> IntervalCenters]

Introduction

Consider the following game. Two players A(be) and B(ert) throw coins. When heads turn up Abe pays Bert a euro, when tails turn up, Bert pays Abe a euro. Before they begin Abe owns € (Subscript[s, A]) and Bert owns € (Subscript[s, B]). The game ends when one of the players runs out of money.

Natural questions are: with what probability will Abe win? How long (i.e. how many throws) is the game expected to last? We can use a Markov chain set-up to find an answer to these questions. In particular we like to investigate how the answers depend on the two parameters (Subscript[s, A]) and (Subscript[s, B]). 

Denote by (Subscript[p, i])(n) with i = 1, 2 ,..., s the probability that Abe owns € i after n coin tosses. The vector with these components is denoted by p(n). In order to reformulate the problem as a Markov chain it is convenient assume that the coin flipping continues, even after the game has ended. Naturally the winner stays the winner, and the loser the loser.

Note for SCI111, 2004: you can skip the following exercise

Exercise 1 (pen and paper)

a) We have p(n+1) = M p(n) for some probability, or transition matrix M. Write down M for the case s = 3.
In the following exercises we take s = 3.
b) There are two absorbing states, which are these?
c) What are the corresponding eigenvectors of M?
d) Assume that A starts with € 2. What is the corresponding start vector p(n) ?
e) Compute the probability that Abe wins in terms of his starting capital (Subscript[s, A]) (so for (Subscript[s, A]) = 1 and (Subscript[s, A]) = 2).

For general s we shall now make a simulation of the game and study the results. For this simulation one very important Mathematica function is the Random function.

throwcoin[p_] := 
        If[Random[] < p, "Heads", "Tail"];

We have used the If statement in the definition. The general format is

	 If[test,t,f] 

where t,f are statements. When this function is called, the test is performed first. This test is function that has True or False as answer.

For example, try

8 < 7
7 < 8
8 == 7

If, in an If statement, test yields True, statement t is performed, if False, statement f is performed.

The throwcoin expression will play the game for us. Try to understand what this function does! Here is an example of its use:

Do[Print[throwcoin[ .5]], {10}]

Repeat the previous command a number of times.

Example: A ball game

We give one example of how to simulate a probability process. This example mainly serves to help you with programming. The problem we consider is a variant of the occupancy problem. We have n urns. In each trial we place one ball in one of the urns. The process ends when one urn contains k balls. This is already quite a complicated process. We would like to have some idea of how long it lasts.

We start by defining n and k.

n = 10;
k = 7;

When the game starts the length of the game is zero and all the urns are empty.

urns = Table[0, {i, 1, n}]
gamelength = 0;

The main loop of the process goes as follows

\(\(While[
        Max[urns] < k, 
\ \ i = 
          Random[Integer, {1, n}]; 
\ \ urns[\([i]\)] = 
          urns[\([i]\)] + 1; 
\ \ gamelength = 
          gamelength + 1
];\)\)
Print[gamelength]

Repeat the two input cells above a number of times and watch how the answer, which is the length of the game, varies. 

What happens here is that we use a While loop.  The While loop

	 While[test,
		statement 1;
		statement 2;
		...
		statement n
	 ] 

executes the statements until test becomes False.

The test is Max[urns]<k. This tests whether in the list urns (that initially just contained zeroes) all numbers are smaller then k. For instance

Max[{0, 5, 6, 8}] < 7

returns False.

The statements that are executed in the body of the While loop represent one trial. In each trial we take a random integer between 1 and  n. 

	Random[Integer,{1,n}];

The outcome is denoted by i, and the number of balls in the i-th urn is increased by 1.

	urn[[i]]=urn[[i]]+1;
	
Finally, we must not forget to increment our trial count

	gamelength=gamelength+1
	
or shorter, in C-style syntax:

	gamelength++
	
Putting it all together in a Module, we get

\(\(\(playthegame[n_, 
        k_]\)\(:=\)\(Module[{urns, gamelength, val, 
          i}, 
gamelength = 0; 
urns = 
          Table[0, {i, 1, n}]; 
While[
          Max[urns] < k, 
i = 
            Random[Integer, {1, 
                n}]; 
\(urns[\([i]\)]++\); 
\(gamelength++\)
]; 
\ \ \ Return[
          gamelength]
]\)\(		\)\)\)

We are ready to do a hundred experiments with for instance n=3, k=7.

res = Table[playthegame[3, 7], {200}]

We can sort the results:

Sort[res]

We can also view the results with the histogram command:

histogram[res]

Simulating the coin game

In the following exercise, but not in the programs to be written, we assume that(Subscript[s, A])+(Subscript[s, B])=3 .

Exercise 2.
a) Write a Mathematica program whowins[sa_, sb_,p_] according to the following specification:

whowins::usage = "whowins[sa,sb,p] returns 1 if Abe wins the game \
and 0 if Bert wins the game. The assets of players Abe and Bert are \
respectively sa and sb. p measures the fairness of the coin, it is the \
probability that heads turn up.";

Hint: The following Mathematica commands (each discussed in this notebook) are likely to be important ingredients:
Module, While, If, Random, ++, --.

b) Let us assume that the coin is fair, i.e. p=0.5. Check your answer on Exercise 1e  by playing the game 10, 100, 1000 times and analyzing the outcomes.
c) Write a howlong program:

howlong::usage = "howlong[sa,sb,p] returns the number of times the \
coin was tossed before the game ends. The assets of players Abe and Bert are \
respectively sa and sb. p measures the fairness of the coin, it is the \
probability that heads turn up.";

d) Again play the game 1000 times with a fair coin (p=0.5) and produce a histogram of the gamelengths.
e) What is the chance that the game is finished in 1 round? And in exactly 2 rounds? And in exactly n rounds?
f) What is the expected number of rounds?
Hint: use Sum.

The formula for success

Now assume that (Subscript[s, A])  +  (Subscript[s, B ]) = 
      10. 
\(The\ coin\ is\ still\ fair . \ 
        We\ are\ again\ interested\ in\ the\ winning\ chances\ of\ \(\(Abe\)\(.\)\)\)

Exercise 3.
a) Play the game 100 times each for every possible asset Abe can start with.
Hint: playing the game once for every possible starting asset could look like this:

s=10;
ex := Table[whowins[i, s-i, .5], {i, 0, s}]

Now Sum the results.

b) Produce a BarChart of the outcomes. (The x axis contains the numbers 0 to 10, indicating Abe's various starting assets. The y axis shows how many of the 100 games Abe won for each starting asset.)

An example of the BarChart command follows:

BarChart[{{1, 3}, {2, 4}}]

Now we move over to general s. But the coin is still fair.

Exercise 4. (Challenge)
a) Prove that the probability that Abe wins when starting with i euros equals

	(i/s)
	
b) Point out the symmetry in this formula.