複数変数を支配するルール

• 一昨日の記事で、n変数をn+1次元空間の円運動で表現する、という話をした
• たとえばロトカ-ヴォルテラに、変数別の周期性が入っているようなデータがあったときに、多変量フーリエ変換をして、その結果を個々の変数ごとに評価する場合と、２変数のスペクトルの強さを合わせて評価する場合とでどんな「スペクトル採用」差が出るかをやってみる

# 多変量フーリエと逆フーリエ
my.mvfft <- function(X,ord.k=2){
	tmp <- mvfft(X)
	md <- Mod(tmp)
	md.sum <- apply(md^ord.k,1,sum)
	ord <- order(md.sum,decreasing=TRUE)
	ret <- tmp[ord,]
	ret.md <- md[ord,]
	ret.md.sum <- md.sum[ord]
	return(list(val=tmp,ord=ord,cmplx=ret,md=ret.md,md.sum=ret.md.sum))
}

my.mvfft.inv <- function(mvout,k=length(mvout$val[,1])){
	tmp <- matrix(0,length(mvout$val[,1]),length(mvout$val[1,]))
	tmp[mvout$ord[1:k],] <- mvout$val[mvout$ord[1:k],]
	ret <- mvfft(tmp,inverse=TRUE)
	tmp2 <- matrix(0,length(mvout$val[,1]),length(mvout$val[1,]))
	for(i in 1:length(tmp2[1,])){
		tmp.ord <- order(Mod(mvout$val[,i]),decreasing=TRUE)
		tmp2[tmp.ord[1:k],i] <- mvout$val[tmp.ord[1:k],i]
	}
	#tmp2[1:k,] <- mvout$val[1:k,]
	ret.conv <- mvfft(tmp2,inverse=TRUE)
	return(list(ret=ret,ret.conv=ret.conv))
}

# ロトカ-ヴォルテラ
# dxdt=ax^k1-bx^k1y
# dydt=cxy^k2-dy^k2
a <- 0.3796097
b <- 0.8440414
c <- 0.5359939
d <- 0.4014024
k1<-2
k2<-3
#a<-runif(1)
#b<-runif(1)
#c<-runif(1)
#d<-runif(1)
#x<-c(runif(1))
#y<-c(runif(1))
# 固定点(x=d/c,y=a/b)の近傍に初期値をとる
x<-c(d/c*0.9)
y<-c(a/b*0.9)

v<-c(-a/y[1]-b*log(y[1])-(c*x[1]-d*log(x[1])))
Niter<-10000
dt<-100/Niter
for(i in 2:Niter){
	dxdt<-a*x[i-1]^k1-b*x[i-1]^k1*y[i-1] 
	dydt<-c*x[i-1]*y[i-1]^k2-d*y[i-1]^k2 
	dxdt <- dxdt*(1+rnorm(1)*100)
	dydt <- dydt*(1+rnorm(1)*100)
	x<-c(x,x[i-1]+dxdt*dt)
	y<-c(y,y[i-1]+dydt*dt)
	v<-c(v,-a/y[i]-b*log(y[i])-(c*x[i]-d*log(x[i])))
}
par(mfcol=c(1,2))
plot(x,y,type="l")

plot(v,ylim=c(-max(abs(v)),max(abs(v))),type="l")
par(mfcol=c(1,1))
tmp.t <- seq(from=0,to=1,length=length(y))*10*pi
n.independent <- 10
for(i in 1:n.independent){
	x <- x+sin(tmp.t*4*pi*runif(1))*runif(1)*0.1
	y <- y+sin(tmp.t*4*pi*runif(1))*runif(1)*0.1

}
plot(x,y,type="l")

mv.out <- my.mvfft(cbind(x,y),2)

par(mfcol=c(2,3))
for(i in 1:length(x)){
	mv.out.inv <- my.mvfft.inv(mv.out,i)
	 matplot(cbind(x,y),type="l",main="original")
	 plot(cbind(x,y),type="l",main="original")
	 par(ask=FALSE)
	 matplot(Re(mv.out.inv[[1]]),type="l",main="楕円合算")
	 plot(Re(mv.out.inv[[1]]),type="l",main="楕円合算")
	 matplot(Re(mv.out.inv[[2]]),type="l",main="通常逆フーリエ")
	 plot(Re(mv.out.inv[[2]]),type="l",main="通常逆フーリエ")
	 par(ask=TRUE)
}
par(ask=FALSE)
par(mfcol=c(1,1))