次元を下げる

• こちらのシミュレーション
• N変量が２つの保存量を持っている
  • １つの保存量は
  • もう１つの保存量は１点を周囲にぐるりと回る保存量
• １つめの保存量から次元は線形にN-1に減らせる
• それをやってみるのが以下のソース

Nx<-5
Nt<-100

T<-sort(runif(Nt))

X<-matrix(runif(Nx*Nt),Nt,Nx)
X[,1]<-X[,1]*10+3
X[,2]<-X[,2]*3*30

MyStandardize<-function(X){
	meanX<-apply(X,2,mean)
	varX<-apply(X,2,var)
	t((t(X)-meanX)/sqrt(varX))
}


###############
n<-5 # 種の数
us<-sample(1:10,n) # それぞれの種の単位Pあたり個体数
# 捕食・被捕食関係の数Nrは種の数と同じ
Nr<-n
# 捕食・被捕食効率
ks<-runif(Nr)
# Inは被捕食種。１度の出会いで殺される個体数
# Playerは捕食・被捕食関係で参加する種とその個体数
# 今は１捕食個体と１被捕食個体の出会いとする
# Outは出会いの結果生じる種とその個体数
# 今は捕食者
In<-Player<-diag(n)
Out<-matrix(0,Nr,n)
for(i in 1:Nr){
	for(j in 1:n){
		if(j==i+1){
			Player[i,j]<-1
			Out[i,j]<-1
		}
	}
}
Player[Nr,1]<-Out[Nr,1]<-1

# ここから時系列シミュレーション
Niter<-10000
dt<-0.01
m<-matrix(0,Niter,n)
p<-rep(0,Niter)
m[1,]<-runif(n)

p[1]<-0
for(i in 1:n){
	p[1]<-p[1]+us[i]*m[1,i]
}
m[1,]<-m[1,]/p[1]
p[1]<-0
for(i in 1:n){
	p[1]<-p[1]+us[i]*m[1,i]
}

for(i in 2:Niter){
	m[i,]<-m[i-1,]
	for(j in 1:Nr){
		tmp<-ks[j]*dt
		pamountIn<-0
		pamountOut<-0
		for(k in 1:n){
			tmp<-tmp*m[i-1,k]^Player[j,k]
			pamountIn<-pamountIn+us[k]*In[j,k]
			pamountOut<-pamountOut+us[k]*Out[j,k]
		}

		for(k in 1:n){
			m[i,k]<-m[i,k]-In[j,k]*tmp+tmp*pamountIn/pamountOut*us[k]*Out[j,k]/pamountOut*Out[j,k]
			#m[i,k]<-m[i,k]-In[j,k]*tmp+Out[j*k]*tmp*pamountIn/pamountOut
		}
	}

	for(j in 1:n){
		p[i]<-p[i]+us[j]*m[i,j]
	}

}
plot(p,type="l",ylim=c(0,2),main="保存量")
plot(data.frame(m),cex=0.1)

#xlim<-ylim<-c(min(m),max(m))
#plot(data.frame(m),cex=0.1,xlim=xlim,ylim=ylim)



##############


X<-m
Xs<-MyStandardize(X)




M<-Xs

# 固有値分解
svdout<-svd(M)
M2<-svdout$u%*%diag(svdout$d) # 分解後データ行列
par(mfcol=c(1,2))
# 固有値分解前後をimage()プロット
image(M)
image(M2)
df1<-as.data.frame(M);df2<-as.data.frame(M2) # データフレーム化
L<-1:5;par(mfcol=c(1,1))
plot(df1[,L],xlim=c(min(df1),max(df1)),ylim=c(min(df1),max(df1))) # ５軸がなす軸ペアでサンプルをプロット。分離しない
plot(df2[,L],xlim=c(min(df2),max(df2)),ylim=c(min(df2),max(df2))) # 固有値分解後に分離力のあるトップ５軸でのプロットは分離する
vM1<-apply(M,2,var) # 項目ごとの分散(分解前)
vM2<-apply(M2,2,var) # 項目ごとの分散(分解後)
ylim<-c(min(vM1,vM2),max(vM1,vM2))
# 固有値分解前の各項目の分散はどれも同じ程度だが
# 固有値分解には、分散の大きいものと小さいものとのコントラストが大きくなっている
plot(vM1,ylim=ylim,type="b")
par(new=TRUE)


plot(vM2,ylim=ylim,type="b",col="red")

par(mfcol=c(1,3))
boxplot(as.data.frame(X))

boxplot(as.data.frame(Xs))


boxplot(as.data.frame(M2))