# Cover.ma
# Cover times for finite state Markov Chains
# by Terry R. McConnell
#
# usage: maple cover.ma (Load, do the calculations, and quit.)
# From the maple command line: read `cover.ma`;
# Load in necessary packages
with(linalg):
with(combinat):
# The following computes the expected amount of time to visit all states
# of a markov chain from the initial state i if the transition matrix is P.
#
# Usage: tcover(Transition matrix, start state);
tcover:= proc(P,i)
local Nn, IndexSet1, IndexSet2, PowerS, Sm, A, PA, n,m, k, Q,
TempSum,j;
# check that the number of arguments is two, and that
# the arguments are of the right type.
if nargs<>2 or not type(P,matrix) then
ERROR(`wrong number of arguments or parameter type in tcover`)
end if:
# check that the matrix is square and the matrix is not empty
# Nn is size of the square matrix
Nn:=rowdim(P):
if Nn<>coldim(P) or Nn<1 then
ERROR(`not a square matrix, or matrix is empty`)
end if:
# we build the index set {1,2,...,Nn}
IndexSet1:={1}:
for k from 1 by 1 to Nn do
IndexSet1:=IndexSet1 union {k}:
end do:
# check for a probablitiy matrix
for m in IndexSet1 do
Sm:=0:
for j in IndexSet1 do
Sm:=Sm+P[m,j]
end do:
if Sm<>1 then
ERROR(`one of rows does not sum to 1`)
end if:
end do:
# we build up the sum in Sm
Sm:=0:
# PowerS is the set of all subsets of {1,2,...,Nn}
PowerS:=subsets(IndexSet1):
# the first set in PowerS is {} and we're not interested
A:=PowerS[nextvalue]():
# for every other subset we have part of the sum
while not PowerS[finished] do
A:=PowerS[nextvalue]():
if not member(i,A) then
n:=nops(A):
PA:=P:
# zero out the rows and columns from the index set A
for k in A do
PA:=addrow(PA,k,k,-1):
PA:=addcol(PA,k,k,-1):
end do:
Q:=inverse( evalm(&*()+ (-1)*PA) ):
# IndexSet2 is the complement of A
IndexSet2:= IndexSet1 minus A:
TempSum:=0:
for j in IndexSet2 do
TempSum := TempSum + Q[i,j]:
end do:
Sm := Sm + (-1)**(n+1)*TempSum:
end if:
end do:
return Sm:
end:
# We apply this to compute the expected time until all possible sequences of
# 0s and 1s of a given length appear in a stream of random digits. In general,
# if k is the length of the sequences sought, then the sequences of that length
# form a Markov chain with doubly stochastic transition matrix P that can be
# easily written down. Indeed, if the states are ordered lexicographically,
# then the first row begins with two consective 1/2's, and each subsequent
# row is obtained by rotating the previous row two postitions to the right.
# The expected time until all sequences of length k have appeared is then
# k + 2^-k Sum tcover(P,x), where the sum is over all sequences x of
# length k.
#
# Here are the results for sequences of length 1:
A:=array([[1/2,1/2],
[1/2,1/2]]);
1+tcover(A,1);
# Below are the commands necessary to do the calculations for length = 2,3,4
# commented out to prevent the calculations from being done when this
# file is loaded.
# Here are the results for sequences of length 2:
# B:=array([[1/2,1/2,0,0],
# [0,0,1/2,1/2],
# [1/2,1/2,0,0],
# [0,0,1/2,1/2]]);
# 2+add(tcover(B,x),x=1..4)/4;
# Here are the results for sequences of length 3:
# C:=array([[1/2,1/2,0,0,0,0,0,0],
# [0,0,1/2,1/2,0,0,0,0],
# [0,0,0,0,1/2,1/2,0,0],
# [0,0,0,0,0,0,1/2,1/2],
# [1/2,1/2,0,0,0,0,0,0],
# [0,0,1/2,1/2,0,0,0,0],
# [0,0,0,0,1/2,1/2,0,0],
# [0,0,0,0,0,0,1/2,1/2]]);
# 3+add(tcover(C,x),x=1..8)/8;
#
# Note that the sum over states is necessary. The time to obtain all
# states starting from state 1 (all ones)
# is longer:
#
# 3 + tcover(C,1);
#
# For sequences of length 4 the computation takes a very long time:
# D:=array([[1/2,1/2,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
# [0,0,1/2,1/2,0,0,0,0,0,0,0,0,0,0,0,0],
# [0,0,0,0,1/2,1/2,0,0,0,0,0,0,0,0,0,0],
# [0,0,0,0,0,0,1/2,1/2,0,0,0,0,0,0,0,0],
# [0,0,0,0,0,0,0,0,1/2,1/2,0,0,0,0,0,0],
# [0,0,0,0,0,0,0,0,0,0,1/2,1/2,0,0,0,0],
# [0,0,0,0,0,0,0,0,0,0,0,0,1/2,1/2,0,0],
# [0,0,0,0,0,0,0,0,0,0,0,0,0,0,1/2,1/2],
# [1/2,1/2,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
# [0,0,1/2,1/2,0,0,0,0,0,0,0,0,0,0,0,0],
# [0,0,0,0,1/2,1/2,0,0,0,0,0,0,0,0,0,0],
# [0,0,0,0,0,0,1/2,1/2,0,0,0,0,0,0,0,0],
# [0,0,0,0,0,0,0,0,1/2,1/2,0,0,0,0,0,0],
# [0,0,0,0,0,0,0,0,0,0,1/2,1/2,0,0,0,0],
# [0,0,0,0,0,0,0,0,0,0,0,0,1/2,1/2,0,0],
# [0,0,0,0,0,0,0,0,0,0,0,0,0,0,1/2,1/2]]):
#
# 4+add(tcover(D,x),x=1..16)/16;