Simulate finite-state process with RandomChoice
- To: mathgroup at smc.vnet.net
- Subject: [mg83054] Simulate finite-state process with RandomChoice
- From: Mark Fisher <particlefilter at gmail.com>
- Date: Fri, 9 Nov 2007 05:09:39 -0500 (EST)
Hi all,
This message is intended to (1) illustrate how to simulate a finite-
state Markov process, (2) highlight a design flaw in RandomChoice, (3)
suggest a fix, and (4) urge the folks at WRI to enhance RandomChoice
to make it more useful.
A transition matrix is a matrix of non-negative numbers where each row
sums to one. Let Q denote a transision matrix. Then Q[[i,j]] is the
probability, conditional on being in state i, that state j occurs
next.
There is no requirement that all of the elements be positive. In fact,
it is often the case that some are zero. However, if all elements are
positive, then it is simple to simulate the process using
RandomChoice. For example, let
Q = {{.2, .8}, {.3, .7}};
NestList[RandomChoice[Q[[#]] -> {1, 2}]&, 1, 100]
For repeated simulations, one may wish to randomize the initial value.
One can use the invariant distribution (if it exists).
InvariantDistribution[Q_] := #/Plus@@#& @
NullSpace[Transpose[Q] - IdentityMatrix[Length[Q]]][[1]]
lambda = InvariantDistribution[Q];
NestList[RandomChoice[Q[[#]] -> {1, 2}]&,
RandomChoice[lambda -> {1, 2}], 100]
However, if any of the elements of Q is zero, then RandomChoice will
not work because it requires all of the weights to be positive. I
think this is a design flaw and I urge the people at WRI to reconsider
this.
For example,
Q = {{.3, .7, 0}, {.6, .1, .3}, {.2, .2, .6}};
lambda = InvariantDistribution[Q];
NestList[RandomChoice[Q[[#]] -> {1, 2, 3}]&,
RandomChoice[lambda -> {1, 2, 3}], 100]
Consequently, one needs to enhance RandomChoice to allow it to work in
this case. Here's one way:
randomChoice::usage = "randomChoice works the same as the built-in \
function RandomChoice except that it allows for non-positive weights."
randomChoice[w:{__?Positive} -> x_, n___] := RandomChoice[w -> x, n]
randomChoice[w_ -> x_, n___] :=
RandomChoice[Rule @@ Transpose[Pick[Transpose[{w, x}], w, _?
Positive]], n]
randomChoice[x___] := RandomChoice[x]
Now we can simulate using randomChoice:
NestList[randomChoice[Q[[#]] -> {1, 2, 3}]&,
randomChoice[lambda -> {1, 2, 3}], 100]
--Mark