Markov Chains - Part 0 - Quick implementation in R

News

Markov Chains - Part 0 - Quick implementation in R

Okay, so doing a lot of work with Markov chains starts to require a bit of optimization. I have found a little trick to get the markov chains work a little faster in R. I will assume you have read about markov chains and had a look at the cool D3 animations here.

So let’s assume we have a state matrix, S, where each row is each individual item and each column is the current state where:

S[x,1..n] (0,1) and S[x,1..n]=1where x (1..m) and x Zand m = number of rows of state matrix and n = number of columns of state matrix

Then we also need a chain matrix, C, where each row contains the probabilities of moving to the next state. Now I know how users think, probabilities are difficult. So let’s say that each row contains a total number of transitions from one state to another. So C is a matrix with n rows and n columns with n defined as from matrix S where:

C[1...n,1...n] Z and C[1...n,1...n] 0

So with my trust issues, let’s create a vector with the totals per row of C so that:

M[1..n] = C[1..n,]

Then we do a element wise division so that each row is divided by the total for the row giving us the probabilities. This allows users to put in counts which is easier to measure in a business process than probabilities and we take care of the conversion.

D=(CT/M)T

Now we get into some of the tricks, change each element is S where it is 1 change it to the column number.

That means S[1..m,1..n] = if 1 then nif 0 then 0

Then use the value in S to create a matrix where each row is replaced with the probability row from C

So the S[1..m,]=C[(S[1..m,]),]

Now we create a random matrix using your preferred random generator either uniform or normal or anything else.

R[1..m,1..m] = Random (0,1) (or if you want to use normal, go for it)

Then you do a piecewise multiplication in R using (*) and not (%*%)

S = R * S

Now all you do is reduce S to only have the maximum value per row and set all other values to 0.

This can all be done without loops in R so that you use compiled code rather than R scripting to do the work.

Here is the code for one loop:

#prepare C

#define C above here

M=colSums(C) D=t(t(C)/M)

#define S as random, just because

S=data.frame(matrix(runif(81*9),81,9))

M=apply(S, 1, max)

S=((S==M)+0)

#one loop through the chain below

#replace the state with the next probabilities

S=D[colSums(colSums(t((S==1)+0)*c(1:ncol(S)))),]

#create random matrix the same size as S

R=matrix(runif(ncol(S)*nrow(S)),nrow(S),ncol(S))

S=S*R

#get the maximum value per row

M=apply(S, 1, max)

#set each row to have 1 in the maximum value’s location, 0 everywhere else.