Simulating Markov Chains In Mathematica
A Markov chain is a sequence of random variables (states) satisfying the Markov property: the probability of the current state depends only on the state that immediately preceded it. In other words, the past state and the future state are stochastically independent. How can we simulate such chains in Wolfram Mathematica?
We define function <$\text{simulateMC}\colon (p_0,p,n) \mapsto (X_1,X_2,\ldots,X_n)$> which will simulate a chain of length <$n$> with initial distribution <$p_0$> and transition distribution <$p$>. Formally, <$p\colon i \mapsto ((p_{i,j_1},\ldots,p_{i,j_k}), (j_1,\ldots,j_k))$> is a function returning the probabilities <$p_{i,j_1},\ldots,p_{i,j_k}$> of states <$j_1,\ldots,j_k \in S$> that can be reached in one step from state <$i \in S$>, where <$S$> is the state space of the chain.
For example, <$p\colon i \mapsto ((\frac{1}{2},\frac{1}{2}),(i-1,i+1))$> is the transition distribution of a simple symmetric random walk (random walks are special cases of Markov chains).
<$\text{simulateMC}$> is implemented in Mathematica as follows:
simulateMC[p0_, p_, n_] := (
chain = ConstantArray[RandomChoice[p0], n + 1];
For[i = 1, i <= n, i++,
chain[[i + 1]] = RandomChoice[p[chain[[i]]]]
];
Return[chain[[2;;]]];
)
Let’s take a look at some examples.
Simple Symmetric Random Walk
We can simulate the simple symmetric random walk mentioned above with this transition distribution function:
randomWalk[pos_] := {.5, .5} -> {pos - 1, pos + 1}
The following piece of Mathematica code will simulate and plot 20 random walks of length 100, each one starting from a random point from <$\{-10,-9,\ldots,9,10\}$>.
ListPlot[
Table[simulateMC[Range[-10, 10], randomWalk, 100], {i, 1, 20}],
Joined -> True
]
Bounded Symmetric Random Walk
Let’s now set up absorbing barriers in <$-10$> and <$10$>, i.e. whenever the walk reaches either of these states, it stays there forever.
boundedRandomWalk[pos_] := If[pos == -10 || pos == 10, {pos}, {.5, .5} -> {pos - 1, pos + 1}]
This time we start all walks in zero:
ListPlot[
Table[simulateMC[{0}, boundedRandomWalk, 100], {i, 1, 20}],
Joined -> True
]
Asymmetric Random Walk
Let’s change the transition probabilities slightly in favor of going up.
asymmetricRandomWalk[pos_] := {.45, .55} -> {pos - 1, pos + 1}
We again start all walks in zero. Note that even a small, 5% change in the probabilities noticeably changes the general tendency of the trajectories.
ListPlot[
Table[simulateMC[{0}, asymmetricRandomWalk, 100], {i, 1, 20}],
Joined -> True
]
Reluctant Random Walk
One more random walk example: transition probabilities now depend on the current position. The closer we get to 0 or 10, the smaller are the probabilities that we are actually going to reach them. The random walk is thus reluctant to end.
reluctantRandomWalk[pos_] := {pos/10, (10 - pos)/10} -> {pos - 1, pos + 1}
We start all walks in 5. Note that the trajectories cluster around 5 which is the mean of 0 and 10.
ListPlot[
Table[simulateMC[{5}, reluctantRandomWalk, 100], {i, 1, 20}],
Joined -> True
]
A Markov Chain Defined By Its Transition Matrix
Finally, let’s simulate a Markov chain having the following transition matrix:
<$$ \left[ \begin{array}{cccc} 0 & \frac{1}{3} & 0 & \frac{2}{3} \\ \frac{1}{3} & 0 & \frac{2}{3} & 0 \\ 0 & \frac{1}{2} & 0 & \frac{1}{2} \\ \frac{1}{2} & 0 & \frac{1}{2} & 0 \end{array} \right] $$>In Mathematica, we can represent this matrix by the following transition distribution function:
chain[i_] :=
Switch[i,
1, {1/3, 2/3} -> {2, 4},
2, {1/3, 2/3} -> {1, 3},
3, {2, 4},
4, {1, 3}
]
We have to reduce the number of chains and their lengths to get a sensible picture.
ListPlot[
Table[simulateMC[{1, 2, 3, 4}, chain, 20], {i, 1, 10}],
Joined -> True
]
Speak your mind
Allowed HTML tags are a, blockquote, em, code, li, ol, p, pre, strong, ul. Links to other comments in the form “[IV]” or “[4]” are detected automatically.