格子座標を正単体座標へ

• k個の量的尺度があって、サンプルについてk個の値のセットが観察されているとする
• サンプルはk次元空間の点に対応させることができる
• こちらでやっているように、代数統計的に考えると、k個の尺度の２値分割表()は、軸の集合のべき集合に分解して考えることができる
• そのやり方をk次元の量的データに持ち込むとしたらどうするだろうか
• サンプルの座標は自由度kであるのに対して、個の頂点を持つ正単体の自由度はであるから、ルールが無ければ１対１対応させることはできない
• k尺度が独立であるものとし、また、各尺度の値を0-1に標準化することが許されるならば、１対１対応の点をとることができる
• ひとたび、正単体座標に展開できれば、複数のサンプルのデータも正単体座標の位置として表現することができる
• これがうまい方法かどうかは、さておき、k次元空間にサンプルが分布する「雲」のようなものが、あるルールで次元空間の１点に対応づけることができるようだ
• この次元空間の点には、各軸、軸のペア、軸のトリオ…というように軸が作るべき集合ごとに情報を取り出す方法が存在する
• 軸ペアの場合は、k個の尺度のペアワイズな相関の強さに対応する(だろう)
• ここで、少し考えてみることがある
• 表では、k個の各軸の比率のみが問題となり、その期待値が、となれば、分散はとなって、分散をわざわざ求める必要もない
• 他方、k軸の量的な値セットの分布を取るときには、各軸に分布があり、n次モーメントもとれるし、モーメントを取るだけでは意味をなさないような分布もあり得る
• すると、べき集合要素が表す「向き」のそれぞれについて分布を考えるときにも、それぞれの「向き」に関する多次モーメントやら、モーメントで表せない分布に依存する何かしら、やらを問題にすることになって、それを問題にすることで、良いことがあるんだろうか…
• というようなことを考えながら、書いたソース(ただのメモ)
  • 多分、こちらの市松模様分布の検出とか、多様体学習とか(こちら)につながる話だし
  • k次元球上の点の位置で表されるベクトルたちの配置の自由度がであって…でも自由度っていう考え方もあって、というあたりの話(こちら)と、一緒に説明されると予想される話

# log-linear modelの期待値の計算関数はちょっと長いですが

LogLinearExp<-function(Table,r,Facets){
	# すべての要素が1である多次元表を作る
	# 期待値表の素
	E<-array(1,r)
	# 表の番地格納した行列を作る
	address<-which(E==1,arr.ind=TRUE)
	# 集合として入力されるFacetsをベクトルのリストにするための格納庫
	FacetList<-list()
	# 各Facetの周辺度数の格納庫
	margList<-list()
	# 周辺度数の番地の格納庫
	margAddress<-list()
	# すべてのFacetsについてループを回す
	# ループ回数のカウンタ
	cnt<-1
	for(i in Facets){
		# Facetの集合を整数ベクトルにする
		FacetList[[cnt]]<-as.integer(i)
		#print(FacetList[[cnt]])
		# このFacetの周辺度数を計算する
		# Facetの要素の数が1であっても、行列として格納されるようにas.matrix()関数を用いている
		if(length(FacetList[[cnt]])==1){
			margList[[cnt]]<-as.matrix(apply(Table,FacetList[[cnt]],sum))
		}else{
			margList[[cnt]]<-as.array(apply(Table,FacetList[[cnt]],sum))
		}
		
		# 周辺度数行列の番地を格納
		margAddress[[cnt]]<-which(margList[[cnt]]<Inf,arr.ind=TRUE)
		#print(margList[[cnt]])
		#print(margAddress[[cnt]])
		# ループ回数を一つ増やす
		cnt<-cnt+1
	}
	# 期待値表のすべてのセルを順次処理するためのループ
	for(i in 1:length(address[,1])){
		#print(address[i,])
		# すべてのFacetに関する処理のループ
		for(j in 1:length(FacetList)){
			# このFacetが集計している軸について、セルの番地を取り出す
			short.address<-address[i,FacetList[[j]]]
			# 軸数が１のときにも行・列の値が必要なのでそのための処理をする
			if(length(short.address)==1)short.address<-c(short.address,1)
			#print(short.address)
			# Facetの周辺度数のうち、該当するものがどれかを調べる
			diff<-apply((t(margAddress[[j]])-short.address)^2,2,sum)
			#print(diff)
			marg.address<-which(diff==0)
			#print(apply(margAddress[[j]],1,identical,short.address))
			#marg.address<-which(apply(margAddress[[j]],1,identical,short.address))
			#marg.address<-which(identical(margAddress[[j]][,1:length(margAddress[[j]][1,])],short.address))
			#print(marg.address)
			#print(margList[[j]][marg.address])
			# すべてのFacetについて積を取る
			E[i]<-E[i]*margList[[j]][marg.address]
		}
	}
	# 期待値表の総数が、元の表の総数と一致するように補正
	#print(sum(Table))
	E<-E/sum(E)*sum(Table)
	# 返す
	E
}
# すべての周辺度数を見てみたいですから、次のようにしてみました
# 上記の処理を関数にする
# install.packages("sets")
CalcAllMarginals<-function(Table,k){
	s<-as.set(1:k)
	powset<-set_power(s)
	marginals<-list()
	cnt<-1
	for(i in powset){
		elems<-as.integer(i)
		if(length(elems)==0){
			marginals[[cnt]]<-list(elems=c(),marginals=sum(Table))
		}else{
			#print(apply(Table,elems,sum))
			marginals[[cnt]]<-list(elems=elems,marginals=apply(Table,elems,sum))
		}
		cnt<-cnt+1
	}
	marginals
}
# 正単体の頂点座標を作る関数
# ncは頂点数
# この関数の説明は面倒くさいけれど…
# 第一ベクトルは(1,0,0,...)として、
# 第二ベクトルは(x1,x2,0,0,...)として
# と、順次、使用軸数を上げて作る
# そのときに「均等に分ける」ことを繰り返す
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])
}
RandomSphere <-
function (df = 3, r = 1, n = 100) 
{
    rs <- matrix(rnorm(df * n), nrow = n)
    rs/sqrt(apply(rs^2, 1, sum)) * r
}

# 周辺度数制約は尺度の組み合わせで決まる
# Facesはそのような尺度の組み合わせの集合(「尺度集合のべき集合」の部分集合)
# そのようなFacesの包含関係を調べて複体表現にする関数

MakeFacets<-function(Faces){
	# Facesの要素である尺度集合同志の真部分集合関係を調べる
	Subs<-outer(Faces,Faces,FUN="set_is_proper_subset")
	# 対角成分は、要素とそれ自身の関係なので、FALSEを入れる
	diag(Subs)<-FALSE
	# ほかの要素に含まれない要素のみが１であるベクトルを作る
	tmp.facets<-abs(sign(apply(Subs,1,sum))-1)
	# ほかの要素に含まれない要素のみをFaces2リストに格納する
	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)
	# Faces2リストの重複を除くために集合扱いする
	tmp<-as.set(Faces2)
	
}
# 分割表の構造nsと周辺度数制約Facetsから
# 分割表のセル数Rに関するRxR行列 X を作るとともに
# 制約がかかっている行のリストを作る
# 尺度ごとに正単体座標を出し
# 尺度の組み合わせは正単体座標ベクトルの外積
# 出力のmatrixXが行列X
# 出力のlistXは行列Xの行ベクトルをリストにしたもの
# 出力のzerosは制約がかかっている行のベクトル

MakeX2<-function(ns,Facets=NULL){
	# 各尺度の正単体座標を出す
	CV<-list()
	for(i in 1:length(ns)){
		CV[[i]]<-CategoryVector(ns[i])
	}
	# すべてのセルの値が0である分割表を作る
	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)
}

# 変換行列Xが回転行列となるようにmatrixXを改変したもの
# listXは改変されていないことに注意
MakeX3<-function(ns,Facets=NULL){
	ret<-MakeX2(ns,Facets)
	X<-ret$matrixX
	dd<-diag(X%*%t(X))
	X2<-diag(1/sqrt(dd))%*%X
	list(matrixX=X2,listX=ret$listX,zeros=ret$zeros)
}
# 分割表の観測度数Obsと分割表の構造nsと周辺度数制約Facetsとから
# 期待値表を返す(返り値のEtable)
# 部分楕球を減次元する行列Xを返す
# 減次元される(Xをかけて0になる)行を返す
# 行列Xは正単体座標とその外積でできていて、回転座標化していないもの
# 多次元市松模様乱点発生

IchimatsuK<-function(n,k){
	r<-matrix(runif(n*k),n,k)
	loop<-TRUE
	ok<-0
	tmps<-NULL
	while(loop){
		s<-matrix(sample(c(0,-1),n*k*3,replace=TRUE),ncol=k)
		tmps<-rbind(tmps,s[which(apply(s,1,sum)%%2==0),])
		if(length(tmps[,1])>=n){
			loop<-FALSE
		}
	}
	tmps<-tmps[1:n,]
	r*((tmps+0.5)*2)
}


k<-3
Npt<-100000
#Ps<-matrix(runif(Npt*k),Npt,k)
Ps<-matrix(0.5,Npt,k)
for(i in 2:Npt){
	Ps[i,]<-Ps[i-1,]
	Ps[i,]<-Ps[i,]+rnorm(k,mean=0,sd=min(abs(c(1-Ps[i,],Ps[i,])))*0.05)
}


library(rgl)
plot3d(Ps[,1:3],xlim=c(0,1),ylim=c(0,1),zlim=c(0,1))

r<-rep(2,k)
R<-prod(r)

# k次元空間の点が観測されたときには、2^k格子点の配置としてk尺度が独立である
# 仮説に基づいて2^k頻度を取らせることにする

Koshi<-function(x){
	tmp<-rbind(1-x,x)
	ret1<-tmp[,1]
	for(i in 2:length(x)){
		ret1<-outer(ret1,tmp[,i],FUN="*")
	}
	ret1
}
Koshi(Ps[1,])

Q<-rep(0,R)
for(i in 1:Npt){
	tmp<-c(Koshi(Ps[i,]))
	Q<-Q+tmp
}
Q.table<-array(Q,r)

Facets<-list()
XZ<-MakeX3(r,Facets)
XZ$matrixX%*%Q/sum(Q)

plot(XZ$matrixX%*%Q/sum(Q))



CalcAllMarginals(Q.table,k)

k<-3
Ps<-RandomSphere(k,0.25,Npt)+0.5
plot3d(Ps[,1:3],xlim=c(0,1),ylim=c(0,1),zlim=c(0,1))

Q<-rep(0,R)
for(i in 1:Npt){
	tmp<-c(Koshi(Ps[i,]))
	Q<-Q+tmp
}
Q.table<-array(Q,r)
Facets<-list()
XZ<-MakeX3(r,Facets)
XZ$matrixX%*%Q/sum(Q)


CalcAllMarginals(Q.table,k)

#########
n<-10000
k<-4
IchimatsuOut<-IchimatsuK(n,k)
library(rgl)
plot3d(IchimatsuOut[,1:3])

plot(as.data.frame(IchimatsuOut))
n<-10000
k<-5
IchimatsuOut<-IchimatsuK(n,k)
plot3d(IchimatsuOut[,1:3])

plot(as.data.frame(IchimatsuOut))

r<-rep(2,k)
R<-prod(r)
Ps<-IchimatsuOut
Ps<-Ps-min(Ps)
Ps<-Ps/(max(Ps))
plot3d(Ps[,1:3])
Npt<-n
Q<-rep(0,R)
for(i in 1:Npt){
	tmp<-c(Koshi(Ps[i,]))
	Q<-Q+tmp
}
Q.table<-array(Q,r)

Facets<-list()
XZ<-MakeX3(r,Facets)
XZ$matrixX%*%Q/sum(Q)

plot(XZ$matrixX%*%Q/sum(Q))

CalcAllMarginals(Q.table,k)