３軸

• カテゴリカルな尺度が３軸ある
  • X,Y,Zの３つ
  • XとYが関係し、XとZが関係している
  • YとZはXとの関係とは独立にさらに関係しているのかどうか知りたい
• ３軸が作る３次元表が観察されたとする
• X,Yの２軸が作る２次元カテゴリの比率は偏りがあるが、そのままにして
• X,Zの２軸が作る２次元カテゴリの比率も偏りがあるが、そのままにして
• Y,Zの２軸が作る２次元カテゴリの比率は、観察データでは偏っているが、「偏っていない」という帰無仮説を想定している
• じゃ、どんな３次元比率が、最尤なのか…を探索してみよう
• 最尤比率が決まれば、自由度(2)で尤度比検定もできる

library(sets)
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])
}
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)
	
}
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
}
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)
}
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)
}
MakeExpected4<-function(Obs,ns,Facets){
	# MakeX2()関数を用いて行列X、制約行を出す
	out.makeX2<-MakeX2(ns,Facets)
	X<-out.makeX2$matrixX
	Z<-out.makeX2$zeros
	# 表の座標にXをかける
	P<-X%*%c(Obs)
	# 期待値表は、制約行が0になっている(それ以外の行はそのまま)なので
	# そのような座標P2を作る
	P2<-P
	P2[which(Z==1)]<-rep(0,length(which(Z==1)))
	# Xの逆行列(solve(X))を用いて、期待値表のセルの値を計算する
	Obs2<-solve(X)%*%P2
	# 期待値表のセルの値はベクトルなので、それを表の構造(array)に納めて返す
	list(Etable=array(Obs2,ns),X=X,Z=Z)
}
MakeExpected5<-function(Obs,ns,Facets){
	out.makeX2<-MakeX3(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)
}

CalcLogLike<-function(Obs,Etable){
	Ptable<-Etable/sum(Etable)
	sum(Obs*log(Ptable))
}


LogLikeRatioMargOK<-function(T1,T2){
	sum(lgamma(T1+1))-sum(lgamma(T2+1))
}

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


library(MCMCpack)

h <- c(rdirichlet(1,rep(1,8)))

array.h <- array(h,rep(2,3))

kajiX <- MakeX2(c(2,2,2))$matrixX

r1 <- r2 <- seq(from = -0.005, to = 0.005, length =100)

r12 <- expand.grid(r1,r2)
LL <- rep(0,length(r12[,1]))
for(i in 1:length(r12[,1])){
	rand.v <- unlist(c(rep(0,6),c(r12[i,])))
	tmp.h <- t(kajiX) %*% rand.v
	new.h <- tmp.h +h
	if(prod(new.h > 0)){
		LL[i] <- sum(h*log(new.h))
	}
	
}

image(matrix(LL,ncol=length(r1)))