simML=function(P,x0,N){
  K=nrow(P)
  x=numeric(N+1)
  x[1]=x0
  for (i in 2:(N+1)) {
    x[i]=sample(1:K,1, prob=P[x[i-1],])
  }
  x
} ## simulira lanac na S={1,2,...,K} s prijelaznom matricom P, 
## x0 je pocetno stanje, N je zadnji trenutak do kojeg simuliramo

P=read.table("matrica_prijelaza.csv", header = FALSE, sep = ",")
# P=read.table("matrica_prijelaza_unif.csv", header = FALSE, sep = ",")
P=as.matrix(P)
P
N=200
set.seed(3024)
traj=simML(P,1, N)
rev_traj=rev(traj)

par(mfrow=c(2,1), mar=c(5,4,1,1))
plot(0:N, traj, type="l", xlab="n", ylab="X_n")
points(0:N,traj, pch=19, cex=0.5, col="blue")
abline(v = 50, col="red")
abline(v = 200, col="red")

plot(0:N, rev_traj, type="l", xlab="n", ylab="X_{200-n}")
points(0:N,rev_traj, pch=19,cex=0.5, col="blue")
abline(v = 0, col="red")
abline(v = 150, col="red")
par(mfrow=c(1,1))

### na jednoj trajektoriji ne mozemo usporedivati razdiobe procesa
### napravimo vise trajektorija:
m=1000 ## broj trajektorija 
       ## (povecati na npr. m=1000 za ilustraciju konvergencije prema stac. distribuciji)
traj_matrica=matrix(NA, nrow=m, ncol=N+1)
set.seed(2242)
for (i in 1:m) {
  traj_matrica[i,]=simML(P,1, N)
}


par(mfrow=c(2,1), mar=c(5,2,1,1))
matplot(0:N, t(traj_matrica), type="l", ylab="")
traj_matrica_rev=apply(traj_matrica, 1, rev)
matplot(0:N, traj_matrica_rev, type="l", ylab="")
par(mfrow=c(1,1))
### ocito nemaju istu razdiobu (gornji krece uvijek iz 1, a donji priblizno iz stacionarne razdioobe
### koja je najvecu masu daje stanju 10)

### trebamo doci do stacionarne razdiobe da bi se vidjela reverzibilnost
par(mfrow=c(2,1), mar=c(5,2,1,1))
matplot(51:N, t(traj_matrica[,51:N]), type="l", ylim=c(1,10), ylab="")
traj_matrica_rev=apply(traj_matrica, 1, rev)
matplot(51:N, traj_matrica_rev[2:151,], type="l", ylim=c(1,10), ylab="")
par(mfrow=c(1,1))


for (k in seq(1,100, by=5)) {
  hist(traj_matrica[,k], prob=TRUE, breaks=seq(0.5,10.5, by=1), 
       main=paste("Distribucija od Xn za n=",k))
  Sys.sleep(1)
}