2011-12-03

分割表の位置・自由度・座標変換５　分割表の正単体・複体座標表現（２）

分割表外積代数複体代数統計 R

できたようだ…
関数

library(sets)
library(igraph)
# 制約条件Facesから、Facets(Facesの包含関係から、最大化した部分正単体の集合としたもの)を作る
MakeFacets<-function(Faces){
	Subs<-outer(Faces,Faces,FUN="set_is_proper_subset")
	diag(Subs)<-FALSE
	tmp.facets<-abs(sign(apply(Subs,1,sum))-1)
	Faces2<-list()
	cnt<-1
	for(i in 1:length(tmp.facets)){
		if(tmp.facets[i]==1){
			Faces2[[cnt]]<-Faces[[i]]
			cnt<-cnt+1
		}
	}
	#tmp<-set(Faces2)
	tmp<-as.set(Faces2)
	
}


# 正単体の座標を出す
# nc個の頂点、nc-1次元
CategoryVector<-function (nc = 3) 
{
    df <- nc - 1
    d <- df + 1
    diagval <- 1:d
    diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - 
        diagval + 2))
    others <- -diagval/(df - (0:(d - 1)))
    m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
    diag(m) <- diagval
    m[upper.tri(m)] <- 0
    as.matrix(m[, 1:df])
}
# 多次元表のdimension nsから
# 表のセル数x表のセル数の変換行列matrixX,(そのリスト型 listX)を作り
# その周辺同数制約Facetsから
# 変換後座標が固定される(自由度が無い)行番号Zを返す
MakeX2<-function(ns,Facets=NULL){
	CV<-list()
	for(i in 1:length(ns)){
		CV[[i]]<-CategoryVector(ns[i])
	}
	t<-array(1,ns)
	taddress<-which(t>0,arr.ind=TRUE)
	R<-prod(ns)
	zeros<-c()
	CVext<-list()
	for(i in 1:length(ns)){
		tmp<-matrix(0,length(CV[[i]][1,]),R)
		for(j in 1:R){
			tmp[,j]<-CV[[i]][taddress[j,i],]
		}
		CVext[[i]]<-tmp
	}


	s<-as.set(1:length(ns))
	ps<-set_power(s)
	X<-NULL
	Y<-NULL
	cnt<-1

	for(i in ps){
		#print(i)
		z<-1
		for(j in Facets){
			if(i<=j){
				z<-0
			}
			#print(i)
			#print(j)
			#print(z)
		}
		if(set_is_empty(i)){
			#print("null")
			tmpX<-rep(1,R)
		}else{
			tmplist<-c()
			for(j in i){
				tmplist<-c(tmplist,j)
			}
			tmplist<-sort(tmplist)
			N<-length(tmplist)
			tmpX<-NULL
			if(N==1){
				tmpX<-CVext[[tmplist]]
			}else{
				initX<-CVext[[tmplist[1]]]
				vals<-NULL
				for(j in 1:R){
					tmpvals<-initX[,j]
					#print(tmpvals)
					for(k in 2:N){
						tmpvals<-outer(tmpvals,CVext[[tmplist[k]]][,j],FUN="*")
						#print(tmpvals)
					}
					tmpvals<-c(tmpvals)
					vals<-rbind(vals,tmpvals)
				}
				tmpX<-t(vals)
			}
			
			
		}
		if(is.null(dim(tmpX)))tmpX<-matrix(tmpX,nrow=1)
		Y<-rbind(Y,tmpX)
		zeros<-c(zeros,rep(z,length(tmpX[,1])))
		X[[cnt]]<-tmpX
		cnt<-cnt+1
	}
	list(matrixX=Y,listX=X,zeros=zeros)
}

# 観察値表と表のdimension,制約条件Facetsから、期待値表を作る
# plot=TRUEの場合には、制約条件に関連してセルに与えられる値をプロットすることによって、制約条件の「イメージ」を確認する
MakeExpected4<-function(Obs,ns,Facets){
	out.makeX2<-MakeX2(ns,Facets)
	X<-out.makeX2$matrixX
	Z<-out.makeX2$zeros

	P<-X%*%c(Obs)
	P2<-P
	P2[which(Z==1)]<-rep(0,length(which(Z==1)))
	Obs2<-solve(X)%*%P2
	
	list(Etable=array(Obs2,ns),X=X,Z=Z)
}
# 多次元表の制約情報であるFacetsを複体(っぽい)グラフにする
GraphFacets<-function(Nv,Facets,plot=FALSE){
	E<-matrix(0,Nv,Nv)
	for(i in Facets){
		i<-unlist(i)
		tmp<-c()
		for(j in i){
			tmp<-c(tmp,j)
		}
		for(j in 1:(length(tmp)-1)){
			for(k in (j+1):length(tmp)){
				E[tmp[j],tmp[k]]<-E[tmp[k],tmp[j]]<-1
			}
		}
	}
	g<-graph.adjacency(E,mode="undirected")
	if(plot){
		plot(g,layout=layout.circle)
	}
	g
}

関数を使ってみる

# 表の次元・次数の設定
Nv<-5
Rs<-sample(2:4,Nv,replace=TRUE)
Vs<-1:Nv
ns<-Rs
# Facesをランダムに作る
Nf<-3

Faces<-list()
maxN<-3
Ns<-sample(1:maxN,Nf,replace=TRUE)
for(i in 1:Nf){
	tmpVs<-sample(Vs,Ns[i])
	Faces[[i]]<-as.set(tmpVs)
}
# Facetsにする
Facets<-MakeFacets(Faces)
# グラフにしてみる

g<-GraphFacets(Nv,Facets,plot=TRUE)

#Facets<-set(1:Nv)
# 観測表をランダムに作る
Obs<-array(runif(prod(ns)),ns)
Obs<-Obs/sum(Obs)

# 期待値表を作る
Etable<-MakeExpected2(Obs,ns,Facets)
Etable2<-MakeExpected3(Obs,ns,Facets)
Etable-Etable2
# ObsとFacetsとから、期待値表を得よう

v2<-c(Obs-Etable)
out.makeX2<-MakeX2(ns,Facets)
X<-out.makeX2$matrixX
Z<-out.makeX2$zeros

P<-X%*%v2
Pz<-X[which(Z==1),]%*%v2

sum((solve(X)%*%P)^2)
sum((solve(X)[,which(Z==1)]%*%Pz)^2)
sum(v2^2)

たくさん回して大丈夫なことを確認しておく
- 観察値と期待値との差ベクトル、それを座標変換したベクトル、座標変換後のベクトルのうち、固定されていない部分空間のベクトル、のそれぞれのノルムが等しいことを、ランダムに表dimension、制約、観察表を発生して確認する

###############xxx###########
Niter<-20
S1<-S2<-S3<-rep(0,Niter)
lenZ<-S1
lenFacets<-S1
lenR<-S1
ngs<-c()
overlap<-rep(0,Niter)
ngZ<-ngP<-list()
cnt<-1
for(ii in 1:Niter){


# 表の次元・次数の設定
Nv<-5
Rs<-sample(2:4,Nv,replace=TRUE)
Vs<-1:Nv
ns<-Rs
# Facesをランダムに作る
Nf<-3

Faces<-list()
maxN<-3
Ns<-sample(1:maxN,Nf,replace=TRUE)
for(i in 1:Nf){
	tmpVs<-sample(Vs,Ns[i])
	Faces[[i]]<-as.set(tmpVs)
}
# Facetsにする
Facets<-MakeFacets(Faces)
# グラフにしてみる
lenFacets[ii]<-length(Facets)
g<-GraphFacets(Nv,Facets,plot=FALSE)
##############
#Facets<-set(1:Nv)
# 観測表をランダムに作る
Obs<-array(runif(prod(ns)),ns)
Obs<-Obs/sum(Obs)

# 期待値表を作る
#Etable<-MakeExpected2(Obs,ns,Facets,plot=TRUE)
#Etable<-MakeExpected3(Obs,ns,Facets,plot=TRUE)
m.e.out<-MakeExpected4(Obs,ns,Facets)
Etable<-m.e.out$Etable
X<-m.e.out$X
Z<-m.e.out$Z
#Etable-Etable2
# ObsとFacetsとから、期待値表を得よう

v2<-c(Obs-Etable)
#out.makeX2<-MakeX2(ns,Facets)
#X<-out.makeX2$matrixX
#Z<-out.makeX2$zeros
lenZ[ii]<-sum(Z)
lenR[ii]<-prod(ns)
#print(Z)
P<-X%*%v2
Pz<-X[which(Z==1),]%*%v2

S1[ii]<-sum((solve(X)%*%P)^2)
S2[ii]<-sum((solve(X)[,which(Z==1)]%*%Pz)^2)
S3[ii]<-sum(v2^2)
ngZ[[ii]]<-Z
ngP[[ii]]<-P
	tmpFlist<-list()
	tmpcnt<-1
	for(iF in Facets){
		tmpFlist[[tmpcnt]]<-iF
		tmpcnt<-tmpcnt+1
	}
	if(length(tmpFlist)>1){
		for(ff1 in 1:(length(tmpFlist)-1)){
			for(ff2 in (ff1+1):length(tmpFlist)){
				if(!set_is_empty(set_intersection(tmpFlist[[ff1]],tmpFlist[[ff2]])))overlap[ii]<-1
			}
		}
	}

if(S1[ii]-S2[ii]>10^(-6)){
	print("-----")
	print(ns)
	print(Facets)
	print(Z)
		print("-----")
	#plot(Z,P)
	ngs<-c(ngs,ii)
}
}

plot(as.data.frame(cbind(S1,S2,S3)))
#plot(as.data.frame(cbind(S1-S2,lenZ,lenFacets)))

この先はこちらへ続く