複素数を使って周回させる - ryamadaのコンピュータ・数学メモ

状態数がk個あり、それらは量を持ち、時間とともに変化するとする
次時刻の状態の量は、現時刻の量から線形に決まるとする
状態推移はkxk行列 M で表される
今、ある状態から始まって、Mがt回適用されたときに、元の状態に戻ることを、周期的とする
周期的な条件は $x =M^t x$ 、すなわち $I=M^t$ である
$M = V \Sigma V^*$ と特異値分解すると $M^t = V (\Sigma)^t V^*$ であるから、 $\Sigma=\Sigma^t$
$\Sigma$ は対角行列なので、対角成分 $S=(s_i),i=1,2,...,k$ が $s_i^t=1$ を満足することがわかる
このような $S$ の条件は、 $s_i$ が絶対値１の複素数であること
以下のように、この条件のみを満足させて作成したMは、要素が複素数の行列となる
ここで、状態の量も複素数とするが、この世で観測できるのは実数部のみであるとする

library(MCMCpack)

k<-30 # No. stata

# 推移行列の特異値を生成(絶対値が１)
e.value<-sqrt(rdirichlet(k,rep(1,2)))*matrix(sample(c(-1,1),k*2,replace=TRUE),k,2)
e.value.Complex<-complex(real=e.value[,1],imaginary=e.value[,2])

# 推移行列のeigen vectorsを生成(ベクトルのModが1)
e.vector<-sqrt(rdirichlet(k,rep(1,2*k)))*matrix(sample(c(-1,1),k*2*k,replace=TRUE),k,2*k)
e.vector%*%t(e.vector)

e.vector.Real<-c(e.vector[,1:k])
e.vector.Imaginary<-c(e.vector[,(k+1):(2*k)])

e.vector.Complex<-complex(real=e.vector.Real,imaginary=e.vector.Imaginary)
e.vector.Complex<-matrix(e.vector.Complex,k,k,byrow=TRUE)

# 検算
m.e.Re<-Re(e.vector.Complex)
m.e.Im<-Im(e.vector.Complex)

m.e.RI<-rbind(m.e.Re,m.e.Im)

t(m.e.RI)%*%m.e.RI

# 推移行列
m<-e.vector.Complex %*% diag(e.value.Complex) %*% solve(e.vector.Complex)

# 初期値を生成
d<-rdirichlet(1,rep(1,2*k))
d.Real<-c(d[,1:k])
d.Imaginary<-c(d[,(k+1):(2*k)])

d.Complex<-complex(real=d.Real,imaginary=d.Imaginary)

# 時間を経過させる
Niter<-100
dseries<-matrix(0,Niter,k)
dseries[1,]<-d.Complex
for(i in 2:Niter){
	dseries[i,]<-m%*%dseries[i-1,]
}
# 一部の状態量の時系列変化をプロットする
tobedrawn<-3
matplot(Re(dseries[,1:tobedrawn]),type="l")