1

I'm wondering if there is an algorithm to simulate a discrete Markov chain with a specific number of occurrences of state knowing the transition matrix way.

For example, how to simulate in R a Markov chain of length n with p occurrences (p < n) of the state '0' for a transition matrix defined by:

TransitionMatrix<- matrix(c(0.7, 0.3, 0.4, 0.6),byrow=TRUE, nrow=2) 


colnames(TransitionMatrix) <- c('0','1')
row.names(TransitionMatrix) <- c('0','1')
Sven Hohenstein
  • 80,497
  • 17
  • 145
  • 168
Armel
  • 31
  • 1
  • 5
  • @nograpes Assuming that _n=20_, `rnd_occ1<-c(0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1)` and `rnd_occ2<-c(0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0)` are two possible candidate having a such transition matrix. Moreover the state "1" has **9** and **6** occurences in the 1st and 2nd sequences respectively. What I'm looking for it is an algorithm or a library to simulate efficiently such a markov chain sequence with for instance _p=12_ occurences of the state "1". – Armel Jan 28 '14 at 18:46
  • I find your comment confusing. You give two possible state transition sequences. Do you want the probability of 12 occurrences of state "1" given a starting state? Do you want to generate random sequences of 0 and 1, and calculate the probability of each of the sequences? It is unclear, and probably effectively answered in the other question. – nograpes Jan 28 '14 at 19:09
  • @nograpes I want to generate random sequences of n numbers (here n=20) of 0 and 1 having p occurences of 1 (here p=12) and knowing the transition matrix of one state to another one. – Armel Jan 28 '14 at 19:24

1 Answers1

0

So, for a two state matrix like yours, this would work to generate random sequences of ones and zeroes, with n total and p ones. It will also give you the probability:

sim.sequence<-function(n,p) {
  ones<-sort(sample(1:n,p,replace=FALSE))
  x<-rep(0,n)
  x[ones]<-1
  row<-as.character(x[1:(n-1)])
  col<-as.character(x[2:n])
  probs<-TransitionMatrix[cbind(row,col)]
  list(states=x,probs=probs)
}

sim.sequence(20,12)
# $states
# [1] 0 0 1 1 1 0 1 0 1 1 0 1 0 1 1 1 1 0 1 0

# $probs
# [1] 0.7 0.3 0.6 0.6 0.4 0.3 0.4 0.3 0.6 0.4 0.3 0.4 0.3 0.6 0.6 0.6
# [17] 0.4 0.3 0.4

Okay, based on your comment, I think I know what you want now. Randomly generated 0's and 1's, and then the proportion of transitions from each state. TransitionMatrix isn't used at all.

So, this would be do it:

sim.sequence<-function(n,p) {
  ones<-sort(sample(1:n,p,replace=FALSE))
  x<-rep(0,n)
  x[ones]<-1
  row<-as.character(x[1:(n-1)])
  col<-as.character(x[2:n])
  tab<-table(row,col)
  probs<-as.matrix(tab) / rowSums(tab)
  list(states=x,probs=probs)
}
set.seed(1)
sim.sequence(20,12)

# $states
# [1] 1 1 1 1 0 1 0 1 1 0 1 0 0 1 1 1 0 0 1 0
# 
# $probs
#    col
# row       0         1
# 0 0.2857143 0.7142857
# 1 0.5000000 0.5000000

I claim that what you are asking for, in the example provided in the comments, is impossible.

You state sequence must have 20 states, meaning 19 transitions. Of those transitions, there are only four possibilities:

(1,1) - m1
(1,0) - m2
(0,1) - n1
(0,0) - n2

So m1+m2+n1+n2==19 and all four are also integers. To restrain b==0.4 it is necessary that m1/(m1+m2)==0.4. Because m1 and m2 also have to be integers, this means that there are only three possibilities for the pair (m1,m2): (2,3), (4,6), (6,9). (8,12) is impossible, because that would imply more than 19 transitions.

Similarly, to restrain a==0.3, you have only one possibility for (n1,n2): (3,7). Therefore, n1+n2==10. This implies that m1+m2==9. However, the only possibilities for m1+m2 is (5,10,15). That is why there will be no string that ever meets your criteria.

Furthermore, we haven't even got into the next constraint where there must be exactly 12 ones. That would make things even more difficult.

nograpes
  • 18,623
  • 1
  • 44
  • 67
  • This is not the solution since the transition probability matrix of the sequence that you generated is matrix(c(1/7,6/7,1/2,1/2),byrow=TRUE, nrow=2) – Armel Jan 28 '14 at 21:07
  • Thanks nograpes. What I want is to randomly generated 0's and 1's under the constraints that (i) the sequence has a length _n>1_ (here n==20) and visits _p_ times (here p==12) the state "1" (ii) the transition matrix of these sequence is `matrix(c(1-a,a,b,1-b),byrow=TRUE, nrow=2)` with _a,b in (0,1)_ (here a==0.3 and b==0.4). So far in what you proposed **only the first constraint (i) is verified**. – Armel Jan 29 '14 at 16:25
  • @Armel What you are asking for is literally impossible, as I show in the comments above. What is really helpful when you ask a question is to provide the expected output, it helps us understand if we understand what you want. – nograpes Jan 29 '14 at 17:51