- こちらから
- 状態数をNとする
- N状態の存在比率をのように時間の関数で表す
- 状態の占拠率の推移をNxN行列で表すとする
- とする
- であるから
- と置けば
- であり、これは
- と解ける
- ただし、
N <- 5
M <- matrix(runif(N^2),N,N)
M <- t(t(M)/apply(M,2,sum))
print(M)
apply(M,2,sum)
TransitionWithExpMat<-function(M,xinit,t,transition=TRUE){
K<-M
if(transition)K<-M-diag(rep(1,length(M[,1])))
e.out <- eigen(K)
V <- e.out[[2]]
U <- solve(V)
s <- e.out[[1]]
X<-matrix(0,length(t),length(M[1,]))
for(i in 1:length(t)){
X[i,]<-(V%*%diag(exp(t[i] * s))%*%U)%*%xinit
}
list(X=Re(X),s=s,V=V,U=U,K=K)
}
xinit<-c(1,rep(0,N-1))
Tmax<-100
L<-1000
t<-seq(from=0,to=Tmax,length=L)
X<-TransitionWithExpMat(M,xinit,t)
matplot(X$X,type="l")