# Stat 753 - Markov Chain Simulation (part II)
# February 18, 2020


# Dimension (number of states)
N=5 

# P is the Transition matrix - randomly generated
P = matrix(0,N,N)
for(i in 1:N) {
  p = runif(N); 
  P[i,] = p/sum(p)
}
P

# Double check that rows of P sum to 1
rowSums(P)

# p is the initial distribution - randomly generated
p = runif(N); p = p/sum(p)
p

# Initial condition for M particles transitioning among the N states
M=1000
x = rmultinom(1,size=M,prob=p) # random initial condition
#x = rep(round(M/N),N) # uniform initial condition

# Do many (100, 1000, more?) iterations using P
iters = 100
for(i in 1:iters) {
  xtmp=matrix(0,N,N)
  for(j in 1:N) {
    xtmp[,j] = rmultinom(1,x[j],P[j,]) 
  }
  x<-rowSums(xtmp)
  barplot(x,main="Current Iteration", ylim=c(0,M/N*1.5), xlab="States", col="seagreen")
  abline(h=M/N) # if random matrix
  Sys.sleep(0.1)
}


# Calculate the stationary distribution
v = eigen(P)$vectors[,1]
v = v/sum(v)
v



###################### Same as above, but here just for 1 particle ########################
# Dimension (number of states)
N=5 

# P is the Transition matrix - randomly generated
P = matrix(0,N,N)
for(i in 1:N) {
  p = runif(N); 
  P[i,] = p/sum(p)
}
P

# Double check that rows of P sum to 1
rowSums(P)

# p is the initial distribution - randomly generated
p = runif(N); p = p/sum(p)
p

# Initial condition for M particles transitioning among the N states
M=1
x = rmultinom(1,size=M,prob=p) # random initial condition

# Do many (100, 1000, more?) iterations using P
iters = 100
for(i in 1:iters) {
  xtmp=matrix(0,N,N)
  for(j in 1:N) {
    xtmp[,j] = rmultinom(1,x[j],P[j,]) 
  }
  x<-rowSums(xtmp)
  barplot(x,main="Current Iteration", xlab="States", col="seagreen")
  Sys.sleep(0.1)
}


# Calculate the stationary distribution
v = eigen(P)$vectors[,1]
v = v/sum(v)
v