複雑な周期性をシミュレーションする

  • こちらこちらで時系列解析について少しまとめた。周期性に関する解析だった
  • 漸化式的な変化の線形代数処理についてはこちらで触れた。時間進行に伴って、確定的だった
  • 順列と置換に関して線形代数的にこちらに書いた。要素を部分集合に分割し、部分集合内での循環を合わせたものであることが示され、それは、特異値分解によって、複素平面上の回転の集まりであることが示された。さらには、そのような置換〜回転は、適当な回数、繰り返すことにより、『元に戻る』ことも示された
  • さて。
  • 独立した循環の繰り返しは、いつか、元に戻ることができることが分かった。
  • 複数の循環があるときに、「交差点」となる要素があって、その要素が、複数の循環に同じタイミングで分配されてしまうと、いつの間にか、「均さ」れてしまう。均されないためには、循環するときに、複数の循環に属してもよいが、個々の時点においては、ただ一つの循環に属するようにすることが大事で、それを守れば、循環のさせ方は、ごちゃまぜでもよいようだ。
  • それを確かめてみる(上段は、短時間にいくつかの要素、下段は、1要素を長時間)
SeqCircuit<-function(s){
	n<-length(s)
	s0<-1:n
	h<-list(c())
	counter<-1
	checker<-rep(1,n)
	while(sum(checker)>0){
		st<-s0[which(checker==1)][1]
		tmpval<-s[st]
		tmpseq<-c(tmpval)
		checker[tmpval]<-0
		while(tmpval!=st){
			tmpval2<-s[tmpval]
			tmpseq<-c(tmpseq,tmpval2)
			tmpval<-tmpval2
			checker[tmpval]<-0
		}
		h[[counter]]<-tmpseq
		counter<-counter+1
	}
	return(list(Circuits=h,NumCircuit=length(h),LengthCircuit=as.numeric(lapply(h,length))))
}

# 要素の量が変化しながら、総量1で固定されている状況

Nelement<-2 # 要素数
# 要素の初期割合
s<-c(rdirichlet(1,rep(1,Nelement)))
S1<-1:Nelement
# 追跡する時間(世代)数
maxT<-500

# 相互に独立するべき循環のセットの数(個々の循環のセットには、複数の循環があるが、それらは、要素を共有しない)
Num<-5

Mlist<-list() # 複数の循環セットを規定する行列を納めたリスト

for(x in 1:Num){

S1<-sample(S1)

tmpm<-matrix(0,length(s),length(s))
for(i in 1:length(s)){
	tmpm[i,S1[i]]<-1
}

Mlist[[x]]<-tmpm
}

# 循環セットの行列を作用させる順番。適当に

rseq<-sample(1:Num,maxT,replace=TRUE)

# 要素の頻度内訳の格納用行列
Ss<-matrix(0,maxT,Nelement)

# 時間発展させる
for(i in 1:maxT){
	s<-Mlist[[rseq[i]]]%*%s
	Ss[i,]<-s
}

# 1からtobedrawnElement個の要素の割合の時間経過をプロット
# 1からtobedrawnTime時間の経過をプロット
par(mfcol=c(2,1))

tobedrawnElement<-min(5,Nelement)
tobedrawnTime<-50
matplot(Ss[1:tobedrawnTime,1:tobedrawnElement],type="l")

tobedrawnElement<-1
tobedrawnTime<-maxT
matplot(Ss[1:tobedrawnTime,1:tobedrawnElement],type="l")

par(mfcol=c(1,1))