推移行列・指数行列

  • こちらから
  • 状態数をNとする
  • N状態の存在比率をP(t)=(p_1(t),...,p_N(t));\sum_{i=1}^N p_i(t)=1のように時間の関数で表す
  • 状態の占拠率の推移をNxN行列M=(m_{i,j});i,j=1,2,...,Nで表すとする
  • P(t+\Delta) = \Delta t \times  M P(t)とする
  • P(t+dt)-P(t) = dt (M - I)P(t)であるから
  • \frac{d}{dt} P(t) = (M - I) P(t)
  • K=(k_{i,j})=M-Iと置けば
  • \frac{d}{dt}P(t) = K P(t)であり、これは
  • P(t) = e^{tK} P(0)と解ける
  • ただし、e^{tK}=(e^{t\times k_{i,j}})
# 状態数
N <- 5
# 推移行列はNxN行列ですべての成分が正、かつ
# すべての列和が1
M <- matrix(runif(N^2),N,N)
M <- t(t(M)/apply(M,2,sum))
print(M)
# 列和が1なことを確認
apply(M,2,sum)
# 推移行列Mと初期値xinitと算出時刻のベクトルtを取って、状態分布行列を返す
# ただし、transition=FALSEとすれば、dx/dt = M x のように微分係数行列を引数とできる
# 返り値
## X は状態分布行列
## sはKの固有値
## V, U は K=V diag(s) Uとなる行列
## K は微分係数行列
TransitionWithExpMat<-function(M,xinit,t,transition=TRUE){
	# K = M-Iを作成
	# ただし、transition=FALSEのときには、MにK=M-Iを取ることにする
	K<-M
	if(transition)K<-M-diag(rep(1,length(M[,1])))
	# Kの特異値分解 ; K = V diag(s) S
	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")